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Abstract 

Geometrically compliant mooring systems that change their shape to accommodate defor¬ 
mations are common in oceanographic and offshore energy production applications. Be¬ 
cause of the inherent geometric nonlinearities, analyses of such systems typically require 
the use of a sophisticated numerical model. This thesis describes one such model and uses 
that model along with experimental results to develop simpler forms for understanding 
the dynamic response of geometrically compliant moorings. 

The numerical program combines the box method spatial discretization with the 
general ized-a method for temporal integration. Compared to other schemes commonly 
employed for the temporal integration of the cable dynamics equations, including box 
method, trapezoidal rule, backward differences, and Newmark’s method, the generalized-a 
algorithm has the advantages of second-order accuracy, controllable numerical dissipation, 
and improved stability when applied to the nonlinear problem. The numerical program is 
validated using results from laboratory and field experiments. 

Field experiment and numerical results are used to develop a simple model for dynamic 
tension response to vertical motion in geometrically compliant moorings. As part of that 
development, the role of inertia, drag, and stiffness in the tension response are explored. 
For most moorings, the response is dominated by inertial and drag effects. The simple 
model uses just two terms to accurately capture these effects, including the coupling 
between inertia and drag. The separability of the responses to vertical and horizontal 
motions is demonstrated and a preliminary model for the response to horizontal motions 
is presented. 

The interaction of the mooring line with the sea floor in catenary moorings is con¬ 
sidered. Using video and tension data from laboratory experiments, the tension shock 
condition at the touchdown point and its implications are observed for the first time. The 
lateral motion of line along the bottom associated with a shock during unloading may be 
a significant cause of chain wear in the touchdown region. Results from the laboratory 
experiments are also used to demonstrate the suitability of the elastic foundation approach 
to modeling sea floor interaction in numerical programs. 

Thesis Supervisor: Dr. Mark A. Grosenbaugh 
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Chapter 1 


Introduction 

A mooring system is typically understood as any type of cable, chain, rope, or tether 
assembly that connects a floating or subsurface buoyant object (ship, buoy, platform) to an 
anchoring system fixed on the sea floor. The floating object will move with environmental 
forcing, but the mooring system will contain the movements to some area (the watch circle) 
centered about the anchoring system. Any mooring system must provide compliance or 
flexibility to accommodate deformations induced by currents and by forcing with periods 
ranging from hours (tides) to seconds (wind waves) without over-tensioning the system 
components. 

This flexibility is typically achieved either through the use of elastically compliant 
members such as rubber tethers or long lengths of synthetic rope, or through geometrically 
compliant configurations in which the system accommodates deformations by changing 
shape without stretching. The geometrically compliant approach is more common in 
situations where adequate compliance or a combination of strength and compliance cannot 
be provided by taught elastic members. This is the case in extremely shallow water, where 
the lengths of the rope or tethers are so short as to limit their compliance. Geometric 
compliance is also often found in offshore deep water applications where the pipe sections 
can be made relatively flexible in bending (through the use of short lengths of pipe and 
flexible joints), but not in axial stretching. Examples of geometrically compliant mooring 
shapes are shown in figure 1-1. 

The shallow water mooring shown in figure 1-1 (a) illustrates the type of mooring typi¬ 
cally used to moor oceanographic, meteorological, and aids-to-navigation buoys in shallow 
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Figure 1-1: Examples of geometric compliance in mooring and riser systems: (a) shal¬ 
low water buoy mooring, (b) deep water oceanographic mooring, and (c) lazy wave riser 
configuration. 

water (on the order of 100 m) [5,13]. The typical mooring for this sort of application 
consists entirely of lengths of chain, with instruments possibly attached between chain 
segments. We say that this type of system is geometrically compliant because its primary 
mechanism to accommodate the motion of the buoy is to lift and lower chain to and from 
the bottom, thus changing its shape. As long as chain remains on the bottom in the 
steady state configuration, the system is typically more flexible geometrically than it is 
elastically. 

The advantages to this arrangement include very high strength due to the use of chain 
as the primary strength member, and the ability to deploy this configuration in a variety 
of water depths. The primary disadvantage to this type of mooring is the need for regular 
replacement of the chain near the bottom of the mooring due to the abrasion of the chain 
on the sea bed [5,23]. A recent alternative to this type of mooring uses elastic tethers as 
the primary compliance mechanism [57,77]. These systems feature significantly reduced 
tensions in most sea conditions because of the much lower mass of the tethers compared 
to chain moorings. Drawbacks to elastic tether moorings include the inability to place 
instruments along the tether and their susceptibility to cutting, either in an accident or 
through vandalism. 

The second type of system in figure 1-1 is an increasingly popular configuration for deep 
water surface moorings for oceanographic applications [31]. Variations on this shape are 
also used to moor meteorological buoys [13]. The s-curve in the mooring shape is achieved 
through careful placement of flotation and ballast along the line. The location of this curve 
at mid-depth allows for geometric compliance without a rigid bottom. Previous deep water 
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surface moorings achieved compliance through the incorporation of long lengths of highly 
stretchable synthetic rope [5]. The advantage to the geometrically compliant system in 
figure 1-1 (b) is the ability to run conductive electromechanical cable along the full length 
of the mooring to bring signals from subsurface instruments to the surface for telemetry. 
Elastically compliant electromechanical members have only recently been introduced [75, 
76] in the oceanographic community and are difficult to handle and relatively expensive, 
particularly for very long lengths. 

Both of the above described mechanisms, an s-shape at mid-depth, and a catenary 
shape along the bottom are often employed together in offshore energy production systems, 
as pictured in figure l-l(c). In this case the mooring line of interest is typically a pipe 
running from the platform to the wellhead. The platform may also be anchored (anchoring 
lines not shown in figure l-l(c)) at multiple points by taut synthetic lines or heavy chain 
and wire lines forming a catenary similar to that described for the shallow water buoy 
mooring. The need for geometric compliance in the riser pipe arises from the inflexibility 
of these pipes to axial (elastic) deformation. 

As a final consideration in this brief overview of compliant moorings, it should be noted 
that in addition to achieving compliance through the mooring line, either geometrically 
or elastically, it is possible in some cases to introduce compliance at the surface by using 
buoys or platforms which have a very low natural frequency (very far below typical wave 
frequencies) such as a spar buoy. This effectively puts a very soft compliant element 
between the wave forcing and the mooring line. Because spar buoys are very long and 
slender, they typically have low reserve buoyancy and are difficult to handle. 

All of the systems pictured in figure 1-1 provide significant compliance to surface wave 
motions under most conditions. One well known mode under which these configurations 
do not provide good compliance is in the case of large currents that pull the geometric 
shaping out of the mooring. The impact of this failure mode can be lessened with the 
addition of elastic compliance into the design. The geometric compliance in these systems 
can also break down during a large storm in which the ability of the mooring to change 
shape may be limited by fluid drag on the cable. This second failure mode is more difficult 
to design for as it can occur even in conditions under which the static shape is preserved 
and may not be alleviated by secondary elastic compliance. Finally, for cases with cable 
resting on the sea floor, friction, adhesion and the elasticity of the bottom can affect the 
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ability of the system to deform geometrically. A loss of geometric compliance as a result 
of any of these mechanisms can lead to dangerously high tensions. A detailed analysis 
of geometrically compliant systems, which will lead to better prediction of these types of 
failures, is the primary goal of this thesis. 


1.1 Analysis of compliant systems 

Much of the recent analytical work relating to geometrically compliant systems has been 
conducted in the context of calculating the contribution of the mooring line damping to 
the overall system dynamics. Brown et al. [7] provide a review of much of the literature 
to date in this area. Most of the work has focused on frequency domain quasi-linearized 
numerical solutions for the slow drift case. Extensive model scale tests have also been 
carried out [55,59,67]. 

Large floating structures (ships, offshore platforms) typically have little damping and 
low natural frequencies for motions in the horizontal plane. For these large structures, 
mooring tensions at wave frequencies are much smaller than the excitation forces. At lower 
frequencies, the mooring forces and wave forces are more comparable. Thus, the damping 
provided by the motion of the mooring system plays a critical role in the response of these 
structures to slow drift motions [59,96]. 

Motions and dynamic tension at wave frequencies are often ignored in these studies. 
This allows for a simplified treatment of the dynamics. For example, Nakamura et al. [67] 
used catenary formulae to calculate the integrated quasi-static velocity and acceleration 
along the mooring. These integrated motions allowed them to write the dynamic tension 
due to slow drift motions in a very simple form. 

In the analyses of compliant systems developed in this thesis, the quantity of interest is 
typically dynamic tension rather than platform motion. Such an approach is particularly 
relevant in oceanographic applications where the motion of the surface platform may not 
be critical, but dynamic tension is dominated by wave induced motions. Knowledge of 
the tension is critical in these applications because components are typically not specified 
with large safety factors for fatigue and ultimate failure (both for cost and ease of handling 
reasons). 

Several authors have considered the impact of wave frequency dynamics on the slow 
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drift damping problem. Huse and Matsumoto [52-55] used a linearized finite element 
model to compute the mooring line damping in the presence of a slow drift regular motion 
superposed with a spectrum of high frequency first-order wave motions. Their calculations 
showed that the damping was two to four times higher when the high frequency motions 
were taken into account. Similar results were obtained by Dercksen et al. [22] and Fylling et 
al. [32] with more sophisticated numerical models. 

In other work that looked at both slow drift and wave frequency excitation, Web¬ 
ster [99] characterized the mooring line damping of a non-dimensionalized catenary riser 
system (a system shaped like that shown in figure l-l(a)) as a function of static tension, 
excitation frequency and amplitude, scope, stiffness, drag, and current. The excitation 
was sinusoidal and either purely vertical or purely horizontal. The numerical model that 
he used was a time-domain nonlinear finite element code. 

Webster [99] also briefly touches on the “impedance” of mooring systems which he 
describes in terms of the trade-offs between geometric and elastic compliance. This is a 
concept first introduced by Triantafyllou et al. [94] to characterize the ratio of elastic to 
catenary stiffness. They noted that fluid drag limits the ability of the mooring to deform 
geometrically and as a result, dynamic tensions increase. 

1.2 Bottom interaction 

An important part of the response in many geometrically compliant systems is the inter¬ 
action of grounded line with the sea floor. Several recent papers have described numerical 
methods for modeling this interaction [16,56,63,89,90]. To date, however, these models 
have not been used to extensively analyze the implications of the bottom interaction on 
the total mooring response. 

Thomas and Hearn [91] and Liu and Bergdahl [63] examined the bottom interaction 
problem in the context of mooring line damping. The results from both papers suggest 
that bottom interaction effects do contribute to mooring line damping, with the in-plane 
friction being more important than the out-of-plane effects [91]. 

Aranha et al. [2], Pesce, Aranha, and Martins [79], and Pesce et al. [80] have examined 
the curvature of riser pipes in the touchdown region using an analytical boundary layer 
approximation. Their goal is to provide better predictions of the bending moment to re- 
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duce fatigue failures. Some of the background for their analytical approach comes from 
work by Burridge et al. [12] and Burridge and Keller [11] for the motion of a string on 
a unilateral constraint. That work demonstrated that a shock wave will form when the 
velocity of the touchdown point exceeds the transverse wave speed of the cable. The ana¬ 
lytical development in Aranha et al. [2] assumes that the touchdown point speed is always 
below this critical limit. No work has been performed that examines the implications for 
mooring dynamics when this assumption does not hold and shock waves do form. 

1.3 Modeling tools 

The problem of predicting the steady state configurations and transient motions of pipe, 
hose, cable, chain, and rope systems in a marine environment is encountered in numerous 
applications. Oftentimes, the methods of solving the problem seem equally numerous. 
Buoy and ship moorings, offshore platforms, and towed systems are often analyzed in very 
different ways, yet are at heart very similar types of structural systems. 

In a 1970 survey paper, Casarella and Parsons [14] compiled an extensive list of work 
related to the hydrodynamic response of cable systems. Their history starts with analytical 
work dating from 1917 to calculate the steady state configuration of cables in air. Through 
1950, treatments of the steady state problem dominated the literature in this area, with 
the first dynamic models for cables in water appearing in 1957. Thomas [90] provides 
a detailed summary of the development of the modern dynamic models, beginning with 
Walton and Polachek’s paper in 1960 [98], and emphasizing developments in the literature 
from the offshore energy field. 

The model developed as part of this thesis provides a nonlinear time-domain solution 
to the mooring dynamics problem. The other modern models described below can be 
similarly classified. Other types of models include frequency domain and linearized or 
quasi-static time domain models. While attractive for their computational efficiency, these 
latter types of models are typically not used for the types of highly nonlinear motions that 
are inherent in the phenomena that are analyzed in the thesis. 
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1.4 State-of-the-art time-domain models for mooring sys¬ 


tems 

Numerical models for mooring systems can be categorized in several different ways. The 
most often cited distinguishing characteristic of a model is the method used to discretize 
the physical system in space. Among the most common methods are finite elements, 
finite differences, and lumped parameter. While there is more universal agreement on 
the temporal discretization method (most use finite differences), there is some variation 
in the way that the temporally discretized equations are integrated in time. Beyond 
these distinctions are the mathematical and physical features incorporated by the various 
models such as bending stiffness, sea bed interaction effects, and treatment of vortex- 
induced vibrations. 

1.4.1 Spatial discretization 

Walton and Polachek [98] published the first treatment of the dynamic solution that re¬ 
sembles very closely the solution methods in use today. They formulate the equations of 
motion for discrete elements and use centered finite differences to discretize the time deriva¬ 
tive terms and step the solution forward in time. With the addition of cable extensibility 
by Polachek et al. [81], a remarkably complete treatment of the nonlinear time domain 
problem existed as early as 1963. This first solution, using a force balance on discrete 
elements to write the equations of motion is what we now categorize as a lumped parame¬ 
ter method. The terminology arises from the lumping of the mass and externally applied 
forces at adjacent nodes which are joined by massless springs. This discretization approach 
has an intuitive simplicity to it and as such is relatively easy to implement. Recent models 
that make use of this approach are described by Huang [47] and Thomas [90,91]. 

In contrast to the summation of forces approach used by lumped parameter methods, 
finite element methods derive their governing equations through principles of virtual work. 
One advantage of this approach is the possibility of a more sophisticated treatment of 
mass. Lumped parameter derivations must necessarily place all mass at discrete nodes 
and then write the governing equations. Finite element methods can derive the governing 
equations using an integration of the mass over the entire element, thus leading to the 
“consistent” mass formulation [62]. The starting point for finite element methods as 
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applied to the marine cable problem is typically a discrete element, much like the lumped 
parameter methods. Examples of such derivations include Engseth [28] and McNamara et 
al. [64], Paulling and Webster [78], following Garrett [33], take the alternative approach 
of formulating differential equations of motion which are solved by the substitution of a 
discrete collection of shape functions which minimize the element energy. The majority 
of state-of-the-art programs currently being used for riser modeling are based on finite 
elements [61]. 

A third approach is to write the continuous partial differential equations and then ap¬ 
ply a spatial discretization scheme based on finite differences. This is the approach taken 
by Ablow and Schechter [1] among others. We distinguish between this and lumped pa¬ 
rameter methods based on the starting point, which in this case is an infinitesimally small 
differential element and in the lumped parameter case is a finite discrete element. Given 
similar physical assumptions the two methods are entirely equivalent, as demonstrated by 
Huang [47]. The distinction between this and the lumped parameter approach is based 
largely on the applications of the method. Many of the numerical solutions for tow cable 
dynamics have used finite differences of the continuous partial differential governing equa¬ 
tions. Another reason for the distinction in this case is simply that, to date, most pure 
lumped parameter methods do not include the effects of bending stiffness in the equations 
of motion [91] 1 . Authors deriving continuous forms of the governing equations have easily 
incorporated this effect [10,46,93]. The model development detailed in chapter 2 is based 
on this approach. 

Finally, a few alternatives to the spatial discretizations outlined above have appeared 
in the literature. Chiou and Leonard [17] and Sun et al. [86] describe the Direct Inte¬ 
gration Method, whereby a boundary value problem is recast as a set of initial value 
problems. Each initial value problem is integrated spatially from a boundary with known 
boundary conditions, and the solutions from these integrations are combined to form a 
total solution that satisfies all boundary conditions. Because the initial value approach 
allows for explicit numerical integration in space, the method has the advantage that the 
solution of large linear systems of equations typical in implicit finite difference and finite 
element schemes can be avoided. There is of course a spatial discretization implied by the 

1 Buckham and Nahon [9] have recently incorporated bending effects into a lumped parameter model for 
low tension ROV tethers. 
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numerical integration of the transformed governing equations. Sun et al. [86] point out 
the need for a method to suppress any spurious solution components that may grow as 
the spatial integrations proceed along the cable. Another alternative scheme is collocation 
which breaks the cable into a small number of segments and fits high order Chebyshev 
polynomials as a solution to the governing equations over each region [15]. 

1.4.2 Temporal discretization 

For all spatial discretization methods the resulting equations are typically written as a 
non-linear matrix equation known as the semi-discrete equation of motion, because the 
time derivatives of the vector of dependent variables are left as continuous functions. 
The exception to this procedure is in finite difference based solutions which typically 
are differenced both in space and in time as part of the same process. This leads to 
yet another distinction between lumped parameter and finite difference approaches. The 
starting point for a finite difference method is typically a set of first-order hyperbolic 
partial differential equations. The equations of motion for lumped parameter schemes 
are most often presented in matrix form as a system of second-order ordinary differential 
equations - the semi-discrete equation of motion. 

Most temporal integration schemes in use today have their roots in the method devel¬ 
oped by Newmark [70], Hughes and Belytschko [50] provide a summary of the development 
of these types of methods in the context of linear finite element structural dynamics. The 
methods typically employ temporal finite differences, with a variety of different schemes 
used to interpolate the solution over the time step. Most classical methods can now 
be cast into unified multi-parameter integration schemes where an adjustment in the 
parameters leads to one of several different methods with different numerical properties 
(e.g., [44,100,102]). Thomas [90] studied the three “classic” methods (Newmark, Houbolt, 
and Wilson-0) and their applicability to the mooring dynamics problem. He concluded 
that Houbolt was the best choice. This is not a surprising result - earlier, Park [74] noted 
that Houbolt was a good choice for highly nonlinear problems. Thomas did not consider 
any of the more modern developments in time integration that are taken up in more detail 
in chapter 2. 

In addition to Newmark and its variants which are popularly employed with finite 
element based models, researchers in the cable dynamics field have employed a variety of 
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different schemes for the temporal integration problem. Chiou and Leonard [17] use simple 
backward finite differences. Sun et al. [86] use the generalized trapezoidal rule which is a 
first-order variant of the Newmark method; it will be discussed in some detail in chapter 2. 
Garrett [33] and Paulling and Webster [78] use the Adams-Moulton method, which in 
first-order form reduces to the trapezoidal rule. Sanders [84] used a computationally 
expensive but fourth-order accurate Runge-Kutta procedure. This is unusual in that most 
researchers have accepted first- or second-order accuracy in order to reduce computational 
expense. 

The most popular finite difference scheme is the box method, in which the governing 
equations are discretized on the half-grid point in both space and time. This method was 
first employed for the solution of tow cable dynamics by Ablow and Schechter [1], Since 
then it has been employed in both towing and mooring applications by Milinazzo et al. [65], 
Howell [46], Tjavaras [93], and Chatjigeorgiou and Mavrakos [15] among others. As will 
be shown, the temporal portion of this discretization is a special form of the generalized-a 
method to be developed in chapter 2. That development will also demonstrate that the box 
method is seldom the best choice of temporal discretization schemes for the cable dynamics 
problem. In a recent paper, Koh et al. [60] came to this same conclusion and proposed a 
modified box method that used backward differences for the temporal discretization. 

1.4.3 Forcing, boundary, and material effects 

There is little disagreement in the proper method of incorporating fluid forces, including 
buoyancy, viscous drag, and added mass forces, into state-of-the-art numerical codes. 
As late as 1970, Casarella and Parsons [14] did choose to distinguish between models 
according to the treatment of drag and whether or not tangential drag was included, but 
there do not appear to be any significant differences between modern approaches. Likewise, 
Breslin [6] laid the groundwork for a consistent treatment of buoyancy and effective tension 
in modern codes. One significant source of hydrodynamic forcing that has not yet been 
fully incorporated into a nonlinear time domain code is vortex-induced vibrations. This 
is an area of active research [95]. 

The numerical treatment of the interaction of the cable with the sea bed is also an 
area of active research. Three basic approaches are prevalent in the literature. Frequency 
domain models (e.g., [94]) and some time domain models (e.g., [89]) cut the mooring off 
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at the touchdown point and attach an equivalent linear spring and/or dashpot. This 
approach is only valid for small dynamic motions about the static touchdown point. A 
second method is the lift-off and grounding approach described by Nakajima et al. [66] 
and Thomas [90]. In this method, the mass of the discrete nodes or elements is reduced to 
zero as they approach the bottom. This simulates a perfectly rigid bottom with no impact 
loads (a smooth rolling and unrolling of the cable, similar to the analytical calculations 
of Aranha et al. [2]). Thomas noted significant numerical difficulties associated with the 
implementation of lift-off and grounding. The third approach is to model the sea bed 
as an elastic foundation. This method has been used by Inoue and Surendran [56] and 
Webster [99]. It is relatively easy to implement and places few restrictions on the types of 
systems that can be modeled. The primary difficulty with this method is in determining 
appropriate elastic and damping constants to associate with a given type of soil. The 
elastic foundation approach is the basis for the bottom interaction model developed as 
part of this thesis. 

For material effects, modern codes may or may not include the effects of material non- 
linearities or bending stiffness. There is little disagreement, however, on the conditions 
under which these effects should be included if an accurate response calculation is to be 
made. Most finite element codes, developed for riser systems that are built from relatively 
large diameter metal pipes, do include bending stiffness, but may neglect material non- 
linearities without a significant loss of accuracy. In the oceanographic community where 
small diameter synthetic mooring lines are common, material nonlinearities can be impor¬ 
tant and bending stiffness can often be neglected. Some codes employ a hybrid approach 
whereby bending stiffness is included only in low tension regions as a numerical smooth¬ 
ing effect (e.g., [87]). A general purpose code should allow for both linear and nonlinear 
materials and for materials with and without bending stiffness. 

1.5 Overview of the thesis 

Chapter 2 describes the development of the generalized-a method for the time integration 
of the cable equations. As an example, the governing continuous partial differential equa¬ 
tions for mooring lines in two dimensions are presented and the reduction to semi-discrete 
form, using spatial finite differences, is derived. The analysis of the stability of a time 
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integration scheme is introduced using the stability of the box method as an example. 
Potential problems with the box method are described and alternative methods are ex¬ 
plored. The generalized-c* method is introduced and the stability and accuracy of the 
method as applied to the cable equations are presented. Comparison is made between the 
new method and many of the previously used methods, including backward differences 
and the generalized trapezoidal family. 

Additional details about the numerical program, including boundary conditions and 
the handling of bottom interaction effects, are described in chapter 3. The algorithms 
used for spatial mesh refinement and adaptive time stepping are also described. Details 
not provided in chapter 3 are given in the appendices. 

The field experiment is described in chapter 4. The centerpiece of the experiment was 
a heavily instrumented all chain mooring. Mooring hardware and instrumentation are 
described. Calibration and data quality issues are also discussed. 

The model is validated and the numerical parameters used in the model are studied in 
chapter 5. The validation is based on analytical and experimental results for a laboratory 
scale hanging chain problem and on full-scale mooring data from the experiments described 
in chapter 4. 

Chapter 6 details a statistical and analytical study of the different contributions to 
the dynamic tension in geometrically compliant systems. These contributions are char¬ 
acterized as drag, stiffness (geometric and elastic), and inertia. The study is based on 
experimental data and extensive numerical runs. Statistical and spectral analyses are 
used along with parametric numerical studies to isolate each of the different tension mech¬ 
anisms. The result of these analyses is a very simple model that can be used to predict 
dynamic tension given a basic characterization of mooring properties, steady state tension, 
and sea state parameters. The chapter concludes with an investigation of the effect of the 
directionality (vertical, horizontal, vertical and horizontal, fully three-dimensional) of the 
input motion. 

A detailed examination of the interaction between the mooring line and the bottom 
is presented in chapter 7. This includes numerical and laboratory simulations of cases 
where there is significant buckling of the line in the region near the touchdown point. The 
implications of the shock condition at the touchdown point are also considered. 

Conclusions and recommendations for follow-on study are presented in chapter 8. 
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1.6 Original contributions 


The numerical program developed in this thesis is based on that of Tjavaras [93] and 
Howell [46]. In the thesis, the program is extended to include bottom interaction effects 
and adaptive discretizations in time step and mesh density. A new temporal integration 
scheme, the generalized-a method, is developed and placed in the context of the recent 
structural analysis literature. An analysis of the stability and accuracy of the overall 
procedure is presented and comparisons are made with other schemes. The new procedure 
has substantially improved stability properties when compared to the old method. The 
model validation detailed in chapter 5 is new for this particular numerical model. 

The analysis of dynamic tension in geometrically compliant systems in chapter 6, 
using regular and random, vertical, horizontal and three-dimensional input motions, and 
a broad range of hydrodynamic and material parameters, is more extensive than any of the 
previous work in this area. The approach to and the development of the simple formula 
for predicting dynamic tension in these system is unique to this thesis. 

Finally, the consideration of the extreme responses of the mooring line on the bottom 
is new. Previous authors [2,94] have limited their analyses to the subsonic case. This is the 
first time that the shock criterion has been experimentally verified and the implications 
of the tension shocks observed and discussed. 
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Chapter 2 


Development of the Time 
Integration Algorithm 


2.1 Governing partial differential equation 

Detailed derivations of the three-dimensional dynamic governing equations for a cable with 
bending stiffness suspended in water are provided by Tjavaras [93]. For completeness, a 
derivation of the two-dimensional equations, upon which the analyses presented in this 
chapter are based, is provided in Appendix A. While the procedure developed below can 
be applied equally well to both two- and three-dimensional models (as will be illustrated 
through the use of both in subsequent chapters), the two-dimensional equations are used 
here for simplicity and succinctness; the two-dimensional model requires six equations 
where the three-dimensional model requires thirteen. The two-dimensional equations for 
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a nonlinearly elastic cable with bending stiffness in steady current are 
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The cable properties are defined by the tension strain relationship, T(e), wet weight, wo, 
mass, m, and added mass, m a , per unit length, diameter, d, and normal and tangential 
drag coefficients, C<j n and Cd t - The motion and force state of the cable is completely 
described by five degrees of freedom (DOF): tangential and normal velocity, u and v, 
strain, e, shear force, S n , and inclination, 4>- A sixth DOF, the curvature of the cable, O 3 , 
is introduced to remove higher order derivatives of 4>■ The current is given in the global 
vertical and horizontal coordinates by U and V, respectively. The relative velocities in 
local coordinates are given by 


u r = u — U cos 4> — V sin 0, 

(2.7) 

v r = v + U sin 4> — V cos (p. 

( 2 . 8 ) 


The independent variables are s, the Lagrangian coordinate measuring length along the 
unstretched cable and t, time. Equations 2.1 through 2.6 can be cast in matrix - vector 
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form as 


3Y dY 

M -^ + K— + F(Y,s,t) = 0 (2.9) 

where Y = [e,S n ,u,v,</),Q 3 ] T and the mass and stiffness matrices, M and K, and the 
forcing vector F are defined in appendix A. 


2.2 Discretization of the governing equation 


The discretization of the partial differential 
governing equation can proceed in several differ¬ 
ent ways. A straightforward method is to use fi¬ 
nite differences in both space and time using the 
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Figure 2-1: Stencil of the box method. 

ences, but because the discretization is applied on 

the half-grid points the method is second-order accurate (see appendix B). The stencil for 
the method is shown in figure 2-1. The result is a four point average centered around the 
half-grid point. Equation 2.9 becomes 


Figure 2-1: Stencil of the box method. 
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The subscripts j define the spatial grid points (the nodes) and the superscripts i define 
the temporal grid points (the time steps). For n nodal points, equation 2.10 defines a 
system of 6 (n — 1) equations to be solved for the 6n dependent variables at time step i. 
The six equations needed to complete the problem are provided by boundary conditions. 
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2.3 Stability of the box method 


The most convenient way to analyze the stability of the box method is to consider the 
stability of the method as applied to an equivalent linear, single DOF system in semi¬ 
discrete form. The first step is to apply the half-grid spatial discretization of the box 
method to equation 2.9. At each half-grid point we derive a set of six equations which we 
can write as 
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where the overdot signifies differentiation with respect to time. The nodal matrices M 
and K, and vector F are defined by 
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0 0 0 


0 

-(m + m a )j -1 
0 
0 
0 
0 


mu + ^ 


a, 1 # 


(nw)j-1 

+ m a ) (U cos <j> + V sin <f>) 

‘ -I 7 — 1 


0 

-(1 + e)j -1 
0 
0 


0 

0 

0 

0 

0 

0 


0 0 —rrij 

0 0 0 

-10 0 
0 0 0 

0 0 0 

0 0 0 


0 

—(m + m a )j 
0 
0 
0 
0 


(mv)j 

mu + ^ p w + ma'j (U cos (j)+V sin 4>) 
0 

~(! + e)j 


0 

0 


0 

0 

0 

0 

0 

0 


( 2 . 12 ) 
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-C3 + 2J-!) 
o 
o 
o 
0 
0 


0 0 0 S n j + S n j_ i 

-2 0 0 -(Tj + Tj-x) 

0—2 0 Vj + Vj-x 

0 0—2 — (uj + Uj_i) 

0 0 0 -2 

0 0 0 0 


0 

0 

0 

0 

-2EI 


T j+ T j- ! 

0 

0 

0 

(Snj + S n j- 1 ) 

0 

0 

2 

0 

0 

Tj + Tj—i 

0 

0 

0 

2 

0 

~(Vj + Vj- 1 ) 

0 

0 

0 

0 

2 

Uj + Uj- 1 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

2EI 


(2.13) 



-Wo (cos (pj +COS0y_i) - ^p w TrdCd, ^(«r |«r| \/l + e)j + ( u r |u r | vTTe) J ._ 1 
Wo (sin (pj + sin <pj-\) + \ p w dC dn [(v r |ry| y/1 + e). + (v r |i> r | y/\ + e) j ._ 1 

0 

0 

~^3 j ~ ^3j-l 

S n (1 + e) 3 j . + [Sn (1 + e) 3 j ^ 


(2.14) 


The shapes of the matrices and vectors in equation 2.11 are diagrammed in figure 2-2. 
If we assemble the blocks associated with the n — 1 nodal matrices and vectors (along 
with appropriate boundary conditions) according to the scheme shown in figure 2-3, then 
it is clear that we can write the semi-discrete equation of motion for all of the dependent 
variables at all of the nodes as 


MY + KY + F = 0. (2.15) 

This is similar to the assembly procedure common in finite element analysis [48]. Prom the 
semi-discrete equation of motion, then, we proceed to reduce the system to a single DOF, 
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Figure 2-2: The shape of the matrices and vectors in equation 2.11. 



Figure 2-3: Assembly procedure for the nodal matrices and vectors into global form. With 
N 1 + N 2 = N total boundary conditions the system is square. The procedure for the global 
mass matrix and force vector is similar. 

linear, homogeneous problem to analyze the stability of the numerical time integration 
procedure. In general, the stability of equation 2.15 in full, nonlinear form, cannot be 
studied analytically. The usual practice is to study the same numerical procedure on 
a simplified model equation, and extrapolate stability properties from there [48, 100]. 
Numerical experiments can then be used to verify the analytical result on the full-scale 
problem. 

The equivalent linear, homogeneous, single DOF problem is 

y + uy = 0. (2-16) 

Applying the box method’s temporal discretization yields a second-order accurate approx- 
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imation for y l 


y l + y l 1 +v(y z + y l l ) =0, 


(2.17) 


where 

y i+iji ~ 1=2 ic~£~)- (2 - 18) 

Rearranging equation 2.18 yields the recursion relationships 

yt = 2 ( ?/ " A f - yi ~\ (2.19) 

yi = y (jP + v^+ir 1 . (2.20) 


If we substitute each of the recursion relationships separately into equation 2.17 we can 
write equations for y l and y l in matrix form as 


yi 


2—uiAt 

2-\-uAt 

0 


yi-l 

_y\ 


-4 

-1 


y i ~ 1 _ 


_ 2+ujAt 



( 2 . 21 ) 


The 2x2 matrix on the right hand side of equation 2.21 is the amplification matrix. The 
spectral radius, p, of this matrix, defined as 


P = max(|Ai| ,|A 2 |), 


( 2 . 22 ) 


governs the growth or decay of the solution from one time step to the next [48]. Ai ; 2 are 
the eigenvalues of the amplification matrix. For p < 1, the solution will remain steady or 
decay and is said to be stable. For p > 1, the solution will grow and is said to be unstable. 
For the time integration scheme defined by the box method, 


2 — tuAf 
1 = 2 + wAf’ 

A 2 = —1, 


(2.23) 

(2.24) 


and the spectral radius is unity (and the scheme is stable) for all values of lo and At. 
When there are no conditions on stability, a procedure is called unconditionally stable. 
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An alternative analysis of the stability of the box method, using classical von Neumann 
stability analysis for finite difference methods, is provided in appendix B. 

In spite of the unconditional stability of the box method, however, the scheme has 
significant problems. Because the update equation for y 1 in equation 2.21 is decoupled 
from y i_1 we can simply write 



/2-^AO 
\2 + uAt) y 


(2.25) 


As LoAt goes to infinity this becomes 



i—1 


(2.26) 


This is the phenomenon known as Crank-Nicolson noise [100], whereby the high frequency 
components of the solution oscillate with every time step. In a linear problem, this noise 
can be removed by computing step-to-step averages once the solution is completed. For a 
nonlinear problem, however, the noise can be a source of instability and hence should be 
eliminated as the solution progresses. 

A second, related, problem is that the spectral radius is constant at unity. An ar¬ 
tifact of the spatial discretization process is that at some point the high frequency (or 
equivalently, high spatial wave number) components of the solution are not well resolved 
and the numerical solution is inaccurate. For this reason it is desirous to have numerical 
dissipation in a scheme such that the spectral radius is less than unity for increasing values 
of u)At [48]. The box method has no numerical dissipation. 

Finally, Wood [100] cites difficulties with averaging schemes in general as applied to 
nonlinear problems. For the nonlinear single DOF case, equation 2.17 can be written as 

f + 1 + ta;V + u/'-y- 1 = 0. (2.27) 


The update equation for y z , equation 2.25, becomes 



V 2 + u'At J y 


(2.28) 


and the stability now becomes conditional as the parameter u changes with time. The 
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practice suggested by Hughes [49], Wood [100] and others, for avoiding this problem is to 
use an averaged value of such as 

■i + -i-l + ^ + ^ 1) = 0. (2.29) 

2.4 Alternatives to the box method 

Given the stability problems associated with the box method, a new solution method 
is sought. Hughes [48] cites the following desirable characteristics in a time-stepping 
algorithm: 

1. Unconditional stability when applied to linear problems: Unconditional stability 
allows the time step to be chosen based on accuracy and resolution concerns, without 
regard for purely numerical issues. 

2. No more than one set of implicit equations to be solved at each time step: This 
minimizes computational expense compared to schemes which may achieve a high 
order of accuracy at a significant computational cost. 

3. At least second-order accuracy: This is a reflection of the constraints imposed by 
Dahlquist’s theorem which states that a third-order accurate method with the most 
appropriate stability conditions does not exist [48]. Again, without a significant 
increase in computational effort, second-order accuracy is the best we can do. 

4. Controllable algorithmic dissipation in the higher modes: In some cases with suffi¬ 
ciently small temporal and spatial discretizations, it may be desirable to have less 
high frequency numerical dissipation. 

5. Self-starting, no information is needed prior to time step zero: Accuracy at time 
step zero (and thus accurate algorithm starting information) is critical in transient 
analysis applications. It is less important in cases where we can slowly ramp up a 
forcing scenario and are not concerned with start-up transients. 

Hulbert [51] adds the following two desirable characteristics: 

6. Single step, that is the solution at i depends only on information at i and i — 1: 
The advantage to a single step algorithm is that it facilitates the implementation 
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of an adaptive time-stepping scheme. If the time step is to be adjusted to A<i in 
going from step i — 1 to step i, then the storage and computational requirement are 
significantly reduced if the solution at i does not also depend on information at i — 2 
which is Ato behind i — 1. 

7. Asymptotically annihilating, or p —> 0 as u>At —> oo: Asymptotic annihilation is 
particularly beneficial in nonlinear problems where it is desirable to damp out high 
frequency noise in just one time step [19]. If the spectral radius at infinity is greater 
than zero, possibly destabilizing noise sources may take several time steps to decay 
completely. 

Finally, based on the idea that nonlinear coefficients should be averaged as discussed 
above, we add that an algorithm should have: 

8. A clear approach to the averaging of temporal coefficient matrices. 

Of unconditionally stable single step algorithms, Thomas [90] compared three his¬ 
torically popular algorithms, Newmark, Houbolt, and Wilson-0, as applied to mooring 
dynamics problems. His conclusion was that Houbolt was the best choice of the three. 
Other recent authors, however, have noted that Houbolt has an undesirable amount of 
low frequency dissipation [19,48]. Also, while asymptotically annihilating, the numerical 
dissipation cannot be controlled (i.e., it can only be asymptotically annihilating). In work 
similar to that described here, Koh et al. [60] proposed a method that retained the box 
method’s spatial discretization but replaced the temporal discretization with a backward 
difference scheme. This scheme is asymptotically annihilating, but only first-order accu¬ 
rate. Sun et al. [86] employ a generalized trapezoidal rule, which does allow for controllable 
dissipation, but is only first-order accurate when there is dissipation. Zueck [103] uses the 
Newmark method, which is the generalized trapezoidal rule for second-order problems, 
and as such also loses second-order accuracy when numerical dissipation is present. 

In the structural dynamics literature, several different schemes have been proposed to 
satisfy the above outlined criteria. Most are developed to solve the second-order semi¬ 
discrete structural dynamics equations, but can be adapted to the first-order problem 
considered here. In fact, equation 2.15 has the same form as the semi-discrete equation 
for transient heat conduction finite element problems. 
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Some of the more recently proposed schemes include the HHT-a [42] and WBZ-a [101] 
methods which combine Newmark style difference formulas with some temporal averaging 
of the terms in the semi-discrete equation of motion. HHT-a averages stiffness, damping 
and force terms. WBZ-a averages the mass terms. Cornwell and Malkus [20] have applied 
the HHT-a method to the first-order semi-discrete heat conduction equation. Bazzi and 
Anderheggen [3] proposed a method whereby the spectral radius at infinity was directly 
set as the sole parameter of the scheme and no coefficient averaging was required. With 
dissipation, however, it is only O(At) accurate. Several multi-parameter “unified” sets of 
algorithms have been published (e.g., [44,71,102]). Through appropriate choices in the 
parameters, these authors are able to implement many of the older methods in addition 
to new schemes. Hoff and Pahl [44, 45] developed what appears to be the most all- 
encompassing such scheme, based on six different parameters. Niemi [71] developed a set 
intended directly for the first-order problem. For our purposes, however, the large multi¬ 
parameter families in their most general forms do not offer a clear and direct approach 
to the temporal averaging of the nonlinear coefficient matrices. A reasonably complete 
family of algorithms that does offer such a clear approach is the generalized-a method 
proposed by Chung and Hulbert [18]. The method is a subset of Hoff and Pahl’s [44] six 
parameter family and can be seen as a straightforward combination of the WBZ-a and 
HHT-a algorithms. 


2.5 The Generalized-ct method 

Cornwell and Malkus [20] applied the HHT-a algorithm to the first-order problem. In this 
method the semi-discrete equation of motion is discretized with temporal averaging of the 
stiffness and force terms, 

MY* + (1 - a)KY i + aKY* -1 + (1 - a)P + aP" 1 = 0. (2.30) 


The difference equation is the same as for the generalized trapezoidal rule [48], 


Y^ 1 + At 


(i- 7 )r- 1 + 7 Y i ' 


(2.31) 
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algorithm 

7 

a k 


1 st order problem 

2 nd order problem 

box method 

1 

2 

1 

2 

1 

2 

Ablow and Schechter [ 1 ] 


backward differences 

1 

0 

0 

Koh et al. [60] 


generalized trapezoidal 

[|,i] 

0 

0 

Sun et al. [ 86 ] 

Newmark [70] 

Cornwell and Malkus 

| — a 

a 

0 

Cornwell and Malkus [20] 

HHT-a [42] 

WBZ-a 

2 + a 

0 

a 


WBZ-a [101] 


Table 2.1: Algorithms included in the generalized-a method. The box method and a 
methods are second-order accurate given the 7 values as shown. The generalized trape¬ 
zoidal rule is second-order accurate only for 7 = |. Backward differences are always 
first-order accurate. 


Following Chung and Hulbert’s development of the generalized-a method for second-order 
equations, we add temporal averaging of the mass terms and equation 2.30 becomes 

(1 - a m )Mr + + (1 - a k )KY i + a^KY ’- 1 + (1 - a k )P + a k P~ x = 0. 

(2.32) 

The three parameter family of algorithms defined by equations 2.31 and 2.32 defines the 
generalized-a method for the first-order semi-discrete problem. Several of the algorithms 
that can be implemented through appropriate choices for 7 , a k , and a m , are summarized 
in table 2.5. 

2.5.1 Accuracy 

As before we analyze the accuracy and stability of the method by studying a single DOF 
problem 

(1 - (x m )y* + a m y I_1 + (1 - a k )u)y z + a k uy l ~ x + (1 - a k )f l + a k f~ x = 0, (2.33) 

y l = y'- 1 + At [(1 - 7 ) 1 + 7 y 1 '] . (2.34) 

The order of accuracy of the method is determined based on a multi-step (information at 
possibly more than just i and i — 1 ), single-stage (only y or y appears, but not both) version 
of the recursion relationship defined by equations 2.33 and 2.34. If we write equation 2.33 
at time step i, eliminate y l using equation 2.34, and add the result to equation 2.33 written 
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at i — 1 and multiplied by At(l — 7 ), we find 


[(1 - a m ) + ju;At(l - ajt)] y l 

- [1 - 2 a m - wA t (1 - 7 - a k + 27 a fc )] y l ~ x - [a m - ujAta k (\ - 7 )] y l ~ 2 
+ At 7 (l - a k )f + At (1 - 7 - a* + 2 7 a fc ) f~ x + Ata fc (l - 7 )f ' 2 = 0. (2.35) 

The local truncation error, r, is the error associated with the use of the difference 
equation 2.35 instead of the exact ordinary differential equation 

y(t) + u>y(t) + f(t) = 0. (2.36) 

If V {?) is an exact solution to this ODE at time t l , then the truncation error is defined 
by [48] 

= i B -y (*'“")] » ( 2 -37) 

n=0 

where B n and C n are the coefficients of the y l and /* in equation 2.35. Expanding y and 
/ terms in Taylor series about t l and then eliminating forcing terms using the exact ODE, 
equation 2.36, yields after some algebraic manipulation 

r (f ) = At (I - 7 - a m + a k ) y (t z ) + O (At 2 ). (2.38) 

Thus, the method is second-order accurate if 

7 + a m -a* = I. (2.39) 


2.5.2 Stability 


Following the same procedure as employed with the box method, the generalized-a method 
for first-order problems can be written in amplification matrix form as 


y l 


V- 1 ' 

= A 

f 


r\ 


(2.40) 
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where the amplification matrix is defined as 


(1 - a m ) + 7(1 - a k )u>At 


(1 - a m ) --fa k coAt At (1 — 7 — a m ) 

-u -a m - (1 -7) (1 - a k )ujAt_ 

(2.41) 


The eigenvalues of this matrix are 


2 [jujAt (a k - 1 ) + a m - 1 ] 


2 a m — 1 + (1—7 — afc + 27 a k ) wA t 


±^ 2 At 2 [(7 - l ) 2 + a k {a k + 27 - 2 )] + 2u>At [7 + 2a m - a k - 1] + 1 . (2.42) 


The method will be stable for all values of ujAt provided that 


<* k < h, < 5, 7 > 5 - 


(2.43) 


Chung and Hulbert [18] suggested a procedure to reduce the scheme to a one parameter 
method. Taking the limit as u>At —> 00, the eigenvalues of the amplification matrix become 


x oo _ f a k 7-1 ) 
1>2 la*-!’ 7 J 


(2.44) 


Requiring second-order accuracy according to equation 2.39 yields A°° as a function of a k 


and a m only 


\ 00 _ J Q fc a fc a m 2 l 

1,2 1 a*-l , a fc -a m + i \' 


(2.45) 


Then, by forcing we can determine a k and a m such that the spectral radius at 

infinity takes on a specific value 


A 00 - 1’ 


3A°° +1 
2A°° - 2' 


(2.46) 


This yields a second-order accurate algorithm in which the only parameter is the spectral 
radius at infinity, p°°. 

Spectral radii of some of the algorithms that are included in table 2.5 along with 
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co-At 

Figure 2-4: Spectral radii of the generalized-a family algorithms. 

results for various values of p°° are shown in figure 2-4. Note that taking A°° G [0,1) as 
the basis for the spectral radius results in a different set of algorithms than A°° £ [—1,0]. 
For p°° = 1 the only option is the negative eigenvalue and this results in the box method. 
A non-dissipative algorithm with A°° = +1 cannot be achieved. The asymptotically 
annihilating form of the algorithm is defined by ak = 0 , a m = — and 7 = 1 . 

The addition of averaging of the mass terms and the a m parameter provides the extra 
degree of freedom that we need to control both the accuracy and the stability over the full 
frequency range. Equations 2.39 and 2.42 define a system of three equations and three 
unknowns. Without the third parameter, a m , we would still have three equations but only 
the two unknowns, 7 and a*. The results from Cornwell and Malkus [20] reflect the fact 
that both Ai and A 2 cannot be controlled while maintaining second-order accuracy. This 
leads to the bifurcations in the spectral radii, evident in figure 2-4, and at some point an 
increase in spectral radius with frequency. Their suggested algorithm is = |, 7 = |. 
Without a m , this is the only possible algorithm that drives the bifurcation point to 00 . It 
is the same algorithm that results from setting A 00 = — | in equation 2.46. 
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2.6 Application to the nonlinear problem 

In applying the generalized-a method to the nonlinear problem we must choose the time 
point at which we will evaluate M, K, and F. A natural choice, consistent with the 
practice suggested by Hughes [49] for nonlinear first-order problems and exemplified by 
equation 2.29, is provided by the temporal averaging of terms that is already a part of the 
method. At time step i equation 2.32 becomes 

(1 - a m )M i ~ am Y i + a m M i_a " 1 Y*~ 1 

+ (1 - ak)K i ~ 0Ck Y i + a k K i ~ 0Ck Y i ~ l + (1 - a k )F + a k F~ l = 0, (2.47) 


where the averaged coefficient matrices are defined as 

M*~ am = (1 - a m )M* + a m M i_1 , and (2.48) 

= (1 - a k )lC + ajfcK'- 1 . (2.49) 

For use with the nonlinear solver described in appendix C, in which the global stiffness and 
mass matrices are never explicitly assembled, it is more convenient to expand equation 2.47 
as 


(1 - a m ) 2 M*Y* + a m ( 1 - a m ) [m^” 1 + M*~ 1 Y i j + a^M^Y^ 1 

+ (1 - a fc ) 2 K i Y i + a fc (l - a k ) [kV ” 1 + K i - 1 Y < ] + a^K^Y^ 1 

+ (1 - a k )F + a k F~ l = 0. (2.50) 

Equation 2.50 represents the temporally and spatially discretized form of the two- or 
three-dimensional cable dynamics equations. The numerical program that implements this 
discretization is described in chapter 3. In chapter 5, this program is used to examine the 
stability of the generalized-a method as applied to the nonlinear cable dynamics equations, 
with particular emphasis on appropriate choices for a k , a m , and 7 (or A^ 2 ). 


50 



Chapter 3 


Implementation of the Numerical 
Program 

The time integration procedure described in chapter 2 is only one piece of the numerical 
program that was developed as part of this thesis. Other important pieces include the 
boundary conditions that round out the governing equations to form a fully determined 
system of equations and the static solution which serves as the initial condition for the 
dynamic solution. The entire solution procedure is diagrammed in figure 3-1. The more 
interesting blocks are described below. Details of the nonlinear solution procedure are 
presented in appendix C. The shooting method solution, which serves as the initial guess 
for the static solution, is described in appendix D. The calculation of coordinate positions 
is presented in appendix E. Details of the program interface and the procedure for model 
and environment description are given by Gobat et al. [35]. 

3.1 Boundary conditions 

As mentioned in the derivation of the semi-discrete equation of motion in chapter 2, the 
governing equations provide only TV x (n — 1) equations for the TV unknown DOF at each 
of the n nodes. The remaining TV equations that are needed to completely determine the 
solution are provided by boundary conditions. The procedures for specifying the boundary 
conditions for the static and dynamic problems are described separately, below. Note that 
much of the complexity in the specification of the static boundary conditions arises from 
the fact that the coordinate positions of the boundaries are not explicitly included as 
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Figure 3-1: Flowchart of the complete numerical solution procedure. Details are given in 
the text and the appendices. 


dependent variables in the governing equations. For a discussion about the merits of this 
formulation see appendix E. 


3.1.1 Static problem 

For the two-dimensional static problem there are four unknowns at each node (N = 4, see 
appendix A for details). The most common boundary conditions are based on specifying 
zero curvature at both ends and applying a known force at the top end. Zero curvature 
is realistic if the cable is attached top and bottom with a joint, shackle, or pivot that 
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releases the moment at the termination. The applied force at the top end comes from 
environmental and other applied forces (a tensioning winch for example) on the platform 


(buoy, ship, drill rig, etc.). The four additional equations are 

flsi = 0, (3.1) 

n 3n = 0, (3.2) 

F x - [T(e n ) cos <f> n - S nn sin <f> n ] = 0, (3.3) 

F y + [T(e n ) sin <j> n + S nn cos (j> n ) = 0, (3.4) 


where F x and F y are the applied forces at node n in the global i and j directions, respec¬ 
tively. 

In many cases, F x and F y are not known directly. For oceanographic surface moorings 
the interaction between mooring forces and buoy forces are coupled through the buoy draft. 
Thus, F x and F y cannot be known before the problem is solved. For offshore applications, 
the specified boundary condition is often the position of the platform relative to the 
anchor and the forces F x and F y are sought as part of the solution. To accommodate 
these conditions we must iteratively solve the static problem with consecutively better 
guesses at the top forces until the desired conditions are satisfied. 

Solving for buoy draft 

Vertical and horizontal forces applied by a surface buoy to the cable segment under the 
buoy are a function of the buoy draft and the known buoyancy and drag properties of 
the buoy. The solution begins with forces calculated from the draft found as part of the 
initial shooting solution, H g . After solving the full nonlinear equilibrium problem, we 
then calculate the actual draft, H®, for these forces based on the position of the top node. 
The absolute error is 

e° H = H° s -H° g . (3.5) 
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that is some small percentage, p ( /, above or below the initial guess, depending on the 
sign of the error. expresses the confidence interval on the initial shooting solution. 
A value of na = 0.1 is typically conservative and works well. With the actual solution 
now bracketed between H ° and Hg, we proceed to use a linear interpolation root finding 
technique [82] to calculate a final solution. This root finding procedure forms a second, 
outer loop of iterations. At each new guessed draft we must go through a new series of 
iterations within the nonlinear solution procedure to solve the problem. The inner loop 
of iterations calculates the equilibrium position for a given applied static force based on 
the current best guess at the draft. The outer iterations continue until the guessed draft 
coincides (to within some specified tolerance) with the calculated draft. 

Resolving platform position 

For the case where we know the position of the upper platform we can use a similar outer 
loop iteration procedure to change the topside applied force until the top end is brought 
into that position. The adjusted applied force at each outer iteration is calculated from 

pk+i = F k -n P ( X k - X) (3.7) 

where F k and X k are the applied force vector and calculated position of the platform 
at outer iteration k, and X is the desired position of the platform. is a numerical 
“stiffness” factor that can be used to accelerate or slow the procedure. These outer loop 
iterations continue until the calculated platform position is within some specified tolerance 
of the known position. The initial values for the forces are determined from the shooting 
method solution which uses this same procedure to bring the platform to the required 
coordinates. 

A third situation requiring outer iterations arises from the inverse of the platform 
positioning problem. In this case, the tension is specified but the horizontal offset of the 
platform relative to the anchor is unknown. In this case we must iterate on the angle 
<j> at the top node such that the specified tension produces forces F x and F y such that 
the platform is on the surface. Like the solution for buoy draft, we can take advantage 
of the fact that the initial shooting solution for (p at the top node should be reasonably 
accurate. Using that initial solution as the first guess, the final solution can be bracketed 
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with a second guess that is only some small distance away from the initial guess. Once <t> 
is bracketed, it can be computed using either bisection or linear interpolation [82], Again, 
the outer iterations proceed until the calculated vertical position of the platform is within 
some specified tolerance of the vertical coordinate of the surface. 

3.1.2 Dynamic problem 

For the two-dimensional dynamic problem with 6 degrees of freedom per node we need 
to formulate a total of six boundary conditions at the two ends. Like the static problem, 
two equations are provided by releasing moments at the two terminations. At the anchor 
we simply impose no motion by setting both normal and tangential velocities to zero. At 
the top we can impose either time varying forces or velocities in the two global directions. 
Velocities are the more common case, as we are typically interested in the response of the 
system to a specified environmentally induced motion of the top of the mooring. In this 
case, the six boundary equations are 


= 0, (3.8) 

u\ = 0, (3.9) 

v\ = 0, (3.10) 

fisl = 0, (3.11) 

Uj - « cos 4> l n - v n sin 4> l n ) = 0, (3.12) 

V} ~ « sin <f> l n + v n cos 4> l n ) = 0, (3.13) 


where Uj and Vj are the specified velocities at time step i in the global vertical and 
horizontal directions, respectively. 

Velocities are typically specified in one of three ways. The first case is a regular motion 
specified as displacements in the two global directions 

xj = A x sin (< jj x t l + ip x ), (3-14) 

yj = A y sin (u> y t l + ip y ) . (3.15) 
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The velocities for this case are 


Uf = A x u x cos (uj x t l + ip x ) , 

Vf = AyUJy COS (Uyt 1 + lj}y) , 


(3.16) 

(3.17) 


where A x%y , u> x ,y > and ipx,y define the amplitude, period, and relative phase of the displace¬ 
ments in the two directions, respectively. 

Secondly, we may specify a random motion profile for a given sea state by breaking the 
spectrum into a summation of individual frequency components with separate amplitudes 
and random phases [29]. For example, a Bretschneider spectrum, specified with a modal 
frequency, w m , and significant height, H s , 


S{u) = 


125 [4 

4 u; 5 


H;e 


2 *- l - 25(^) 4 


(3.18) 


can be discretized over m frequencies, with a spacing of Aw. The amplitude of the 
k th component is 


A k = V2S{u k )Au>. (3.19) 

The displacement is the sum of all the discrete components 
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A i = ^ A k sin (tf +if>k)- (3-20) 

k=1 

The random phases, ipk, are generated as uniform random numbers on the interval [—7r, 7 t]. 
The total velocity is 

m 

U} = ^2 A k u k cos (u k t % + ip k ). (3.21) 

k=1 

This procedure is not limited to spectra which are known analytically. It can easily be 
applied to wave spectra derived from field data gathered by such instruments as wave¬ 
following buoys and acoustic doppler current profilers. 

Finally, for model validation purposes it is often convenient to impose an entire time 
series of motion onto the top of the mooring. These time series might be the integrated 


56 


motions from accelerometer data that were recorded during storm events. Given the known 
platform motion we can compare model predicted tensions to those actually recorded in 
the field. The velocity record necessary for this application can either be numerically 
integrated from acceleration or numerically differentiated from displacements, depending 
on the available data. If the velocity record consists of discretely sampled points, with 
a spacing between points of A t v then the velocity at time step i is interpolated by 

U} = (U^ 1 - 1 7*) (^- - k + l) + 0* (3.22) 


where k defines the appropriate index into the zero-offset velocity record, 



(3.23) 


3.2 Bottom interaction 


Following the same basic approach as Webster [99], the unilateral boundary condition 
at the sea floor is modeled as an elastic foundation with linear stiffness and damping 
properties. Given the vertical coordinate of the bottom, which may vary with horizontal 
position, Xbottom(?/)) the bottom exerts a force on node j if Xj < bottom (%')• For both 
static and dynamic problems the force is defined as 


F b = k\ Xj \, (3.24) 

where k is the stiffness per unit length of the bottom. In static problems the force is 
constrained so that Fj, < wo- The force is always assumed to act in the global vertical 
direction and as such can be treated simply as a modification to the wet weight, wo, in 
the governing equations. In the dynamic problem we also add a damping force, 


F d = —bvj, (3.25) 

where b is the dashpot constant of the bottom and Vj is the normal velocity of node j. 

One of the disadvantages of this approach is that appropriate values for k and b are 
difficult to calculate without extensive field and laboratory experimental testing of soils. 
For most problems, however, the gross response of the system is largely insensitive to the 
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choice of values. Typically, we specify k as the fraction of the line wet weight that will 
be supported with a deflection equal to one diameter. A non-dimensional form of the 
stiffness, k, can be defined as that fraction, 

k = — - (3.26) 

wo 

The damping constant b is calculated from a specified value of a damping ratio, (. Given 
£, the mass plus added mass of the grounded line, m + m a , and the natural frequency of 
the elastic foundation/cable system, o; n , the damping constant is [92] 

b = 2( (m + m a ) co n . (3.27) 


The natural frequency of the system is calculated as 






k 


m + m a 


(3.28) 


A damping ratio of 0.5 is typically sufficient to eliminate any spurious high frequency 
effects that result from the line impacting the bottom without significantly affecting overall 
system response. 

The advantages to this treatment of the bottom are the simplicity with which it can be 
implemented and the generality which it allows. The approach places no restrictions on 
the number of touchdown points or where and how those points move during the dynamic 
problem because the entire mooring, including grounded line, is always “in play”. This 
contrasts with approaches which may track a single touchdown point, adding or removing 
line from the problem to calculate a dynamic response only for line that is instantaneously 
above that point. The implementation described above has no difficulty handling cases 
in which positively buoyant portions of the line float above the bottom between heavier 
portions of line which remain on the bottom or in which a traveling wave of ungrounded 
line moves along a portion of grounded line. 
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3.3 Refinement of the spatial discretization 


In many moorings with low flexural stiffness, the half grid spatial discretization can lead to 
undesirable spatial oscillations in the solution. This phenomenon can be easily understood 
by considering the equation relating shear force to curvature, 

EI^ + S n ( 1 + e ) 3 = 0. (3.29) 

For a static solution this equation is discretized as 

2 El — J -— + S n j (1 + ej) 3 + S n j_i (1 + ej_i ) 3 = 0 . (3.30) 

/\Sj J 

If El w 0 and e << 1 as is typical, the only solution (barring A Sj = 0) is S n j & — S' n ,-_ 1 . 
If IS'nl > 0, the shear force will oscillate about zero from one node to the next. This error 
is particularly manifest in areas of high curvature and at the boundaries. The problem 
can be minimized by increasing the density of the spatial mesh [10]. 

Without a priori knowledge of the static solution the most easily applied spatial dis¬ 
cretization is uniform, 


Asj — 


L 

n — 1 ’ 


(3.31) 


where L is the length of the cable segment and n is the number of nodes used in the 
discretization. Unfortunately, a uniform mesh with small As to reduce spatial oscillations 
can require large numbers of nodes. An alternative is to make the mesh finer only in 
problem areas: areas of high curvature and at the boundaries. To automate this allocation 
of nodes we can develop a procedure such that given a static solution based on a uniform 
mesh, we can optimize the mesh in some sense and then recalculate the static solution to 
take advantage of the refinement. The procedure outlined below is based on that described 
by Eggleton [27]. It is worth noting that Press et al. [82] describe a procedure, also based 
on Eggleton’s approach, that adaptively refines the mesh as part of the nonlinear solution 
procedure. That procedure had significant problems with convergence when applied to the 
geometrically nonlinear problems considered here. It also requires that three equations 
and additional dependent variables be added into the problem. 
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The approach to mesh refinement can be understood by considering a minimization of 
the sum given by 

n 

Y. [c w {^z{ s j) — ^z( s j-i)) + ( s j ~ s j-i)] 2 • (3.32) 

3 =2 

The Sj coordinates of the n nodes are unknown, but from our previous static solution we 
can provide a good estimate of fl 3 (s) for any value of s. The first term in the sum will keep 
nodes close together in areas of high variability in fl 3 . The second term will keep nodes 
from getting too far apart in areas with low variability in f 2 3 . c w is a constant that controls 
the weighting used to place nodes with respect to the two effects. c w » 1 will result in a 
large proportion of the n available nodes being used in areas of high curvature with large 
spacing between the remaining nodes in other regions of the system. In contrast, c w « 1 
results in a nearly uniform mesh, with little emphasis placed on refining mesh density in 
high curvature regions. 

Minimizing equation 3.32 requires that we solve an n degree of freedom nonlinear least 
squares problem. Alternatively, we can approximate the sum as an integral and cast the 
minimization as a variational problem. If we define the mesh control function 

/(s) = c«,fi 3 (s) + s, (3.33) 


then equation 3.32 is simply 

£[/(<¥)-/(*!- 1 )] 2 - (3-34) 

i=2 

Without affecting the solution of the minimization problem we can introduce a new in¬ 
dependent variable, q, that varies uniformly throughout the mesh (A q is a constant) and 
rewrite the summation as 


E 


3 =2 


fjsj) - /(gj-i) 


2 



s j s j— 1 


A q 


(3.35) 
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We can approximate this sum as the integral given by 



The integral form can now be minimized by solving the variational problem 

6 k,fr s (s) + ll 2 — ds = 0. 

Jo dq 


(3.36) 


(3.37) 


In writing equation 3.37 we have substituted a normalized estimate of the curvature gra¬ 
dient, Sl' s , for the spatial derivative of Q ! s is defined as 



1°3 j ~ ^3j-l 
max | 0 3fc - k _i 


, k — 1 ... n 


- 1 . 


(3.38) 


This formulation normalizes the curvature to have a maximum value of one and a lower 
bound at zero. c w can then be interpreted as the mesh density weight for curvature effects, 
relative to unity. 

The solution to the variational problem in equation 3.37 can be written as [21,27] 


ds _ (3 

dq c w £l' s + 1 ’ 


(3.39) 


where f3 is a constant to be determined. Equation 3.39 is a boundary value problem for 
s with boundary conditions s = 0 at q = 0 and s = L at q = L. We use the shooting 
method [82] and bisection to determine (3 such that all boundary conditions are satisfied. 
Bounds on (3 can be derived by considering the extreme cases fl((s) = 0 and 0((s) = 1; 
both conditions lead to a uniform mesh, 


f 2 ],(s) — 0 —> — — 1 — 1 >• Anin — 1> (3.40) 

ds 

0 (,(s) = 1 — > — = 1 Anax = 1 + Cy,. (3-41) 


With each trial (3, we integrate from q = 0 to q — L using fourth order Runge-Kutta 
integration. The error function for the bisection is simply s(q = L) — L, i.e., the difference 
between the integrated s coordinate of node n and its known coordinate, s n = L. 

The final step in the process it to recalculate the static solution on the refined mesh. 
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With some care, it is possible to minimize the computational expense associated with 
outer iterations during this second solution because the solution for the unknown boundary 
conditions is unlikely to vary significantly from the uniform mesh solution to the refined 
mesh solution. 


3.4 Adaptive time-stepping 

The stability analysis of the generalized-a method that was presented in chapter 2 can be 
strictly applied only to a linear form of the problem. In the nonlinear case the method 
cannot guarantee stability because the nonlinear solution procedure at each time step 
is not unconditionally convergent. Because the nonlinear solver uses the result from the 
previous time step as the initial guess at the solution for the current time step, the solution 
may not converge if those two solutions are significantly different. For this reason, there 
are limits to the maximum allowable time step that can successfully be used to propagate 
the solution in time without giving rise to numerical instabilities. 

Typically, we choose a value of At based on factors such as the accurate resolution 
of the physics in the problem and the desired sampling rate of the numerical solution. 
Depending on the particular problem this value of At may not be small enough to avoid 
numerical instabilities that arise over the course of the simulation. This situation is 
common in cases where the cable goes slack for brief periods of time or when there is 
rapid lifting and lowering of cable to and from the bottom. A procedure for avoiding the 
numerical problems in these cases, without modifying the baseline time step for the whole 
problem, is adaptive time-stepping. 

The adaptive time-stepping procedure that is implemented here is relatively simple. 
If an instability arises the time step will be reduced automatically to try to get through 
that portion of the simulation. At each time step where the baseline time increment is not 
small enough to accurately propagate the solution, At is reduced by a factor of ten. The 
solution then proceeds through ten steps at the smaller increment. The reduction can be 
recursive, with a practical limit set as four orders of magnitude below the base value of At. 
If the nonlinear solver fails even at this lowest increment, the solution is aborted. This 
procedure has the advantage that the simulation always produces results on the originally 
requested sampling grid. 
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Adaptive time-stepping is only of limited usefulness, however, without some care being 
taken in the choice of a baseline time increment. If the algorithm is deciding that it needs 
a smaller time increment at every step then it would be faster to have set a smaller time 
step in the first place. 
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Chapter 4 


Description of the Field 
Experiment 


Numerical studies of complex mechanical systems, like the geometrically compliant moor¬ 
ings considered in this thesis, have the advantages that they place few constraints on the 
system under study and that they are relatively inexpensive to conduct. In contrast, exper¬ 
imental efforts are limited by practical and cost considerations. Nevertheless, a numerical 
study alone is seldom able to paint a complete picture of the physics that govern the re¬ 
sponses of these kinds of systems. For this reason, both a field and a laboratory experiment 
were conducted as part of this thesis to provide an added level of detail and confidence. 
The field experiment described in this chapter offered the chance to collect full scale data 
that reflect a response to real environmental conditions. The results from the experiment 
are used in chapter 5 as part of the validation of the numerical program and in chapter 6, 
along with simulation results, to analyze the dynamic response of mooring systems. The 
laboratory experiments described in chapter 7 provide higher spatio-temporal resolution 
under more controlled conditions. These advantages facilitate the detailed analysis of the 
bottom interaction described in that chapter. 


4.1 Location and climatology 

The Shallow Water Engineering Experiment (SWEX) was conducted at an area known as 
the WHOI Buoy Farm. This is a one km 2 area approximately 40 km southwest of Woods 
Hole, Massachusetts or 18 km southwest of Gay Head on the island of Martha’s Vineyard, 
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Figure 4-1: Geographic map of field experiment site. The star marks the WHOI Buoy 
Farm. The base GPS receiver station was located at the Gay Head lighthouse on the 
southwestern most tip of Martha’s Vineyard, marked by the black square. The Buzzards 
Bay tower is marked by the black circle. 

Massachusetts. The site location within Rhode Island Sound is shown in figure 4-1. The 
locations of the moorings within the Buoy Farm are shown in figure 4-2. Nominal water 
depth at the site is 42 m. 

The experiment was deployed on 5 December 1998 and recovered on 20 January 1999 
to coincide with a portion of the winter storm season. As shown in figure 4-1 the site 
is exposed with significant fetch to wind and storms from the south, southeast, and to a 
more limited degree the southwest. There is much less exposure to significant storms from 
the north and northeast due to limited fetch. Based on climatological records from the 
nearby Buzzards Bay C-MAN tower, the dominant winds blow from the southwest during 
this period. Figure 4-3 shows the hourly averaged wind records from the Buzzards Bay 
tower during the experiment. Through December and January the average wind direction 
was 224° (southwest) and average wind speed was 18.3 knots. There were several large 
storm events, however, with winds from the southeast. The largest of these occurred on 
3 January 1999, with peak sustained winds of 50 knots. 
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Figure 4-2: Location of the experimental moorings within the Buoy Farm during SWEX. 
Surf = surface mooring, ST = Seatex waverider buoy. The 600 kHz ADCP was located in 
a bottom mounted tripod that was on the groundline between SSB and SSB P/U. Dashed 
circles indicate the approximate watch circle of each mooring. A,B,C, and D mark the 
four corner guard buoys that delimit the Buoy Farm. Other markers indicate additional 
experiments and fishing floats that were deployed at the Buoy Farm during the field 
experiments. 

4.2 Mooring hardware 

The primary experimental mooring was an all chain catenary mooring. The mooring 
design is shown in figure 4-4. The system consisted of 80 m of |-inch galvanized steel 
trawler chain, broken only by three inline accelerometer instruments (AxPacks). The 
AxPacks were hose clamped onto stainless steel strongbacks and the strongbacks were 
shackled between shots of the chain. The surface buoy was a cylindrical block of Surlyn 
foam 1.27 m in diameter and 0.75 m high. An instrument well extended through the 
middle of the buoy and 1.4 m beyond the bottom of the foam. The well is approximately 
24 cm in diameter. The properties for all of the mooring components are summarized in 
table 4.1. 
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Figure 4-3: Winds during the field deployment. The data are hourly averaged results 
from the Buzzards Bay C-MAN tower maintained by the National Data Buoy Center 
(http://www.ndbc.noaa.gov). Shaded areas indicate dates where all channels (dark) or 
just the y accelerometer channel (light) of the surface buoy instrument had significant 
data errors. 


material 

length (m) 

m (kg/m) 

w 0 (N/m) 

EA (N) 

Tj-inch chain 


3.73 

31.85 

6.4 x 10 7 

AxPack 

0.76 

10.02 

70.82 

8.0 x 10 7 

shackle/ring/shackle 

0.22 

16.22 

81.23 

8.0 x 10 7 


Table 4.1: Properties of the components used in the experimental mooring. AxPack 
properties include two |-inch chain shackles and a |-inch pear ring at each end of the 
strongback. The axial stiffness of components that include a shackle are based on the 
stiffness of a shackle. 
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= 5/8" chain shackle 
3/4” pear ring 
5/8" chain shackle 


(IT) = 5/8" chain shackle 
3/4" pear ring 
1" anchor shackle 



2300 lb anti-fouling anchor 


Figure 4-4: Schematic of the surface mooring used in the field experiments. 


4.3 Instrumentation 

4.3.1 Engineering instrumentation 
Mooring line instrumentation 

The mooring chain was instrumented with three AxPack self-contained accelerometer in¬ 
struments as shown in figure 4-4. They were located so as to span the region of high 
curvature near the touchdown point over the range of currents that were expected at the 
site. The lowest instrument was placed so that it would be approximately 3 m off the 
bottom at the lowest tide and slack current. The data indicate that it probably did hit 
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the bottom at various times during the deployment. 

Each AxPack consists of Tattletale Model 8 microcontroller (with eight channels of 
12-bit A/D) from Onset instruments mated to a Persistor CF8 compact flash board (with 
a 24 MB compact flash card) from Peripheral Issues. The accelerometer is a Summit 
Instruments model 34103A triaxial accelerometer with a 0 - 5V output scale over the 
range ±1.5(7. The primary advantage to these accelerometers is that they are completely 
self-contained. Given a single +5V power supply they produce an amplified and filtered 
0 - 5V signal. The internal filter is a single pole Butterworth filter with the 3 dB point 
at 4.6 Hz. The accelerometer is packaged in a small cube less than 2.5 cm on a side. 
Power is provided by three 3.6V lithium C cells (manufactured by Tadiran). All power 
conditioning is done on board the Tattletale. 

The sample rate throughout the experiment was 10 Hz. The AxPack accelerometers 
were sampled for 20 minutes beginning on the hour at 0800, 1600, and 0000 hours localtime. 
Because there is no communication between the instruments during the experiment, each 
unit carries a separate battery backed real time clock (Dallas Semiconductor DS1302). 
These clocks were synchronized using an electronic trigger pulse prior to deployment. The 
crystals for these clock chips appear to have been cut from the same batch and exhibit 
similar drift characteristics, with each AxPack losing approximately 50 seconds in 30 days. 
These clocks retain the real time in case of a fault and reset in the Tattletale. 

The electronics and accelerometer are secured into a machined aluminum rack and 
together with the batteries sealed into delrin pressure cases. The pressure cases are 21 cm 
long and 7.5 cm in diameter. A photo of the assembled and unassembled AxPack com¬ 
ponents is shown in figure 4-5. The driving factor in the AxPack design was to keep 
the size, mass, and wet weight of the units as consistent as possible with the rest of the 
mooring. However, on their strongbacks and taking into account the shackle/ring/shackle 
assemblies that are required to attach the AxPack inline with the rest of the mooring, 
the AxPacks have approximately twice the mass and wet weight per length of the ^-inch 
trawler chain. 

Buoy instrumentation 

The buoy was instrumented with a six axis motion package: triaxial linear acceleration 
(Columbia Research Laboratories model SA-307-TX) and three Systron Donner single axis 
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Figure 4-5: AxPack strongback, pressure case, and electronics. 


gyro rate chips. The instruments were controlled and logged (at 12-bit resolution) by a 
Tattletale Model 6F controller. The surface package also included a Precision Navigation 
TCM2 electronic compass module. The digital signal from the compass was converted to 
an analog signal using the onboard 8-bit digital to analog converter. This analog signal 
was then sampled by the Tattletale for logging, providing 256 levels of heading around the 
360° of the compass in the final dataset. The connection to the mooring chain was made 
through a 5000 pound capacity load cell. The load cell was also sampled by the Tattletale. 
All of these instruments were sampled at 12.5 Hz (though the effective update rate of the 
compass is only 1 Hz) three times a day for twenty minutes. Due to a programming 
error prior to deployment, the start time of each sample was delayed by five minutes 
relative to the AxPack sample periods; the three sample periods began at 0805, 1605, and 
0005 localtime. No attempt was made to synchronize the surface instrument clock with 
the AxPacks beyond setting them within approximately one second of each other before 
sealing the instruments. 

All of the instruments performed well for the first three weeks of the deployment. Data 
from the surface buoy instruments had significant drop-outs and obvious signal problems 
from 27 December through 31 December. After 31 December, the y accelerometer signal 
(one of the horizontal axes) was always bad, but the other channels appeared to be problem 
free. During a post-deployment analysis it was determined that the multiplexer channel 
for the y accelerometer had failed. Our speculation is that while it was in the process of 










failing it caused problems with the other channels, but that once it had failed completely, 
the remaining channels were unaffected. 

GPS instrument 

The surface buoy also contained a GPS (global positioning system) receiver. The moti¬ 
vation for including this instrument was to verify the quasi-static position of the buoy 
within its watch circle. From the ship’s GPS during deployment we knew the location 
of the anchor to within several meters. By recording GPS signals at the buoy and at a 
non-moving base station located at the Gay Head lighthouse on the island of Martha’s 
Vineyard (figure 4-1) we hoped to resolve the motions of the buoy to within better than 
one meter [24]. The GPS receivers were Canadian Marconi Allstar units with 1 Hz po¬ 
sition, velocity, and time output and 1 Hz carrier phase output. On the buoy the GPS 
receiver was controlled by a Tattletale Model 8 with logging to a Peripheral Issues Per- 
sistor AT8 with a 175 MB flash card. The base station GPS receiver was controlled by a 
standard PC. Unfortunately, the remote receiver failed. We feel confident, however, that 
the technique can provide an interesting and valuable dataset and thus the system will be 
redeployed on a future engineering test deployment. 

4.3.2 Environmental instrumentation 

In order to quantify the environmental forcing on the surface mooring, both waves and 
current were measured during the deployment. Current was measured using two acoustic 
doppler current profilers (ADCPs): a 600 kHz unit mounted in a tripod on the sea floor 
and a 1200 kHz unit mounted on top of a subsurface buoy that was tethered at 13 m 
depth. Directional wave spectra were measured by a Seatex Wavescan buoy (Seatex A/S). 
This buoy is moored such that it has a significant portion of its tether floating on the 
surface. This allows it to respond freely to the incident waves in heave, pitch, and roll. 
The motion of the buoy is measured using a six axis Hippy unit. 

As part of a separate effort, the Wavescan data will be compared with the ADCP 
data to test the ability of the ADCP to resolve directional wave spectra. This comparison 
required relatively high frequency sampling of the ADCP. The ping rate was 3 Hz, with 
the velocity results averaged and stored at 1 Hz (i.e., each 1 Hz current sample represents 
the average current result from three pings over the previous second). The current data 


72 



is provided as a profile, with 75 cm between depth bins on the 600 kHz unit and 25 cm 
between bins on the 1200 kHz unit. Accounting for the tripod height above the bottom, 
the first current point is 1.95 m above the bottom. An overly conservative number of bins 
was used so that the last bin always fell beyond the surface. The ADCPs were sampled 
for 40 minutes (600 kHZ) and 26 minutes (1200 kHz) three times per day (0800, 1600, 
0000 ). 

4.3.3 Data telemetry 

All of the instruments stored their data locally. The Wavescan buoy and the surface buoy 
both had ARGOS satellite transmitters that were used for location purposes only. This 
allowed remote monitoring of the location of the buoys to ensure that they had not failed 
and gone adrift. 


4.4 Data processing 

All accelerometer and gyro calibrations were performed using manufacturer supplied cali¬ 
bration coefficients. The validity of the accelerometer calibrations was verified both before 
and after the deployment through a check of each instrument’s outputs in a variety of 
tilted positions. The 5000 pound load cell was sent to the manufacturer for a recalibration 
immediately prior to the experiment. 

The motion of the buoy in earth referenced coordinates was computed using the ap¬ 
proach outlined by Edson et al. [25]. In this procedure, the orientation of the local 
coordinate system is computed using a complementary filter in which the high frequency 
signal from the rate gyros is combined with low frequency tilt and heading information 
derived from the horizontal accelerometers and the compass. 1 The result of the comple¬ 
mentary filter is a time series of buoy orientation which can be used to transform the 
recorded accelerometer signals into east, north, and vertical components. These earth ref¬ 
erenced accelerations are then integrated into velocity and displacement, with a highpass 
filter at each step to remove any low frequency (greater than 30 second period) drift. 

1 Results after 27 December 1998 were processed with the assumption that the low frequency y accelerom¬ 
eter signal was identically zero, i.e., that there was no systematic tilt in that direction. The assumption 
is easily justified given the near zero mean y accelerometer signals prior to 27 December and it allowed 
us to compute an estimate of the vertical motion even after the loss of the y accelerometer. Motions in 
the horizontal plane (east and north) were not computed for data after 27 December. 
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For use with the numerical model it is convenient to determine an approximate ori¬ 
entation for the plane of the mooring and to project wind, current, and motion vectors 
into a coordinate system oriented with this plane. This approximation and projection 
allows us to apply forcing data derived from the three-dimensional experimental results in 
two-dimensional model simulations. Definitions for the procedure are shown in figure 4-6. 
We determine the direction of the plane of the mooring by considering the net force due 
to wind and current on the buoy only. We neglect for now any current drag on the chain 
because the currents tended to decay sharply away from the surface and thus drag forces 
on the chain were much smaller than the forces on the buoy due to current and wind. For 
a current profile V(z) with magnitude and direction at the surface Vh and 9h, and wind 
with speed W and direction 8 W , the north and east components of the force, F n and F e , 
are 


F n 


cos 6h cos 6 w 


\pS b c db v 2 

A 


sin 6h sin 0 W 


\p*rS w C dw W 2 _ 


where S w , Cdwi Sb and Cdb are the projected area and drag coefficient above and below the 
buoy waterline, respectively. The resultant effective direction of the plane of the mooring 
is 

6 e g = tan -1 • ( 4 - 2 ) 


Given the effective plane determined by 6 e g, we seek effective values of the wind, W e g, 
and current profiles, V e g(z), which yield the same level of force as the true forces projected 
onto the effective plane 


F c cos (0 e ff ~ 6{z)) = \pS b C db V e g(z ) 2 , (4.3) 

F jj COS {&eff 9w) = 2 Pair S w C d w ff'eff ' ( 4 - 4 ) 
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Figure 4-6: Definitions for the procedure to determine the approximate two-dimensional 
plane of the mooring. 

Taking care of signs, we define the effective current and wind in this plane as 

V eS (z) = sign [cos (0 e ff - 0(z))} V ( z ) \/|cos (0 eff - 0{z))\, (4.5) 

W eff = sign [cos (<9 e ff - 0 W )] W \/|cos (0 eff - 9 W )\. (4.6) 

The north and east components of the buoy motion, X n (t ) and X e (l), respectively, are 
converted into in-plane and out-of-plane components according to 


Afp — X n cos 0 e jj + X e sin0 e ff, (4-7) 

X 0 p — X n sin $eff “h Xe COS $eff • (4.8) 

Because the average water depth of 42 m was near the maximum range of the 600 
kHz unit, the data from the 1200 kHz unit appear to be more accurate near the surface. 
As the near surface currents (along with the wind) dominate the steady-state response of 
the mooring, the profiles, V(z ), used in the procedure outlined above are based on data 
from the 1200 kHz instrument, with extrapolated values below 13 m depth. While both 
instruments recorded three ping ensembles at 1 Hz, only temporally averaged profiles (over 
the 26 minute length of each sample period) were used. 
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Chapter 5 


Validation and Parameterization 
of the Numerical Program 


The numerical model described in chapters 2 and 3 and appendices C through E was 
validated using the data collected during the full scale field experiment and by comparison 
with two hanging chain problems with known solutions. The first step in the validation 
is to characterize those aspects of the model which are purely numerical, particularly the 
time integration and mesh refinement parameters. We do this by comparing simulation 
results to a known solution. This allows us to establish reasonable values for the numerical 
parameters which are then used in the comparison of model results with experimental 
results in order to establish the ability of the numerical program to accurately predict 
dynamic response under a variety of forcing conditions. 


5.1 Response of a hanging chain 

Figure 5-1 depicts the hanging chain system used for the first part of the validation. Two 
cases will be considered. In the first case we apply a very small initial displacement to the 
chain and then at time t = 0, release it. The dynamic response of the chain for t > 0 can 
be calculated analytically for the small motions that result. In the second case we impose 
a sinusoidally varying horizontal displacement to the top of the chain and analyze the 
forced response. This latter problem was studied both numerically and experimentally by 
Howell [46]. 
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5.1.1 Free response to an initial displacement 

For small motions and an inextensible chain, the equation of motion for the chain can be 
written as 


d 2 q d f dq 
= m9S ¥s 


(5.1) 


Assuming a harmonic solution of the form 


q(s,t) = q(s ) [Acoscvt + B sin uit] 


(5.2) 


the mode shapes, q(s), are [97] 


q(s) = ci J 0 



+ C2Y0 



(5.3) 


The requirement that the solution be finite at s = 0 leads to the elimination of the Yo 
term and the requirement that q(L) = 0 leads to the natural frequencies, u>. They are 
given by the roots of 


Jo 




(5.4) 


The complete response is the sum of the response in each mode 


q(s, t) = ^ Jp ^2 oj n [^n cos ut + B n sina;t]. 


(5.5) 


The coefficients A n and B n are determined from the initial displacement, qo(s), and 
velocity, qo{s). Given qo(s) — 0, we can immediately determine that B n — 0. To determine 
A n we first write 


?( 5 >0) = ^ A„J 0 ^2u; n y^ = q 0 (s). 


(5.6) 


Multiplying both sides by Jo ^2o>„ ) > integrating from s = 0 to s = L, and making use 
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of the fact that 



For comparison with simulation results, the analytic 
solution was computed for a chain released from an initial 
catenary configuration. For simplicity all of the model 
parameters (mass per length, gravity, length) were set to 
unity. The horizontal force applied at s = 0 to create 
the initial deflection was set to 0.001. For simulation 
results El was set to 10~ 6 and EA to 10 9 . All of the 
integrals were computed using the trapezoidal rule with 
10000 panels. A 400 second time series of the response Figure 5-1: Definitions for the 
at the free end was constructed using the first 20 modes hanging chains problems, 
of the analytic solution. The analytic result was sampled at 20 Hz; the natural frequency 
for mode 20 is approximately 5 Hz. 

Because the primary distinction amongst the various algorithms contained within the 
the generalized-a method is the amount of numerical damping, all results are compared 
in the frequency domain. For each 400 second time series, power spectra of the response 
at the free end were computed using non-overlapping 256 point FFTs. As an example, 
figure 5-2 shows the power spectra for the analytic solution and for a numerical solution 
with At = 0.01 s, and 50 nodes. The circle and square markers indicate the 

spectral peaks as computed using a crude peak detection algorithm. In subsequent results 
only the peaks are plotted. 

Figure 5-3 shows a comparison between the analytic solution and numerical solu¬ 
tions for six different parameterizations of the generalized-a method. At this time step, 
At = 0.01 s, most of the algorithms are accurate out to the fifth or sixth mode. The 
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Spectral Response Comparison 



Figure 5-2: Power spectral density of the response of the free end of the chain for the 
analytic solution and for a numerical solution with A^ 2 = — §, At = 0.01 s, and n = 50. 

notable exception is the first-order accurate backward differences, which substantially un¬ 
derestimates the response even in the first mode. All of the algorithms show a marked 
fall off from the analytic solution at higher frequencies, with the solutions for A^ 2 > 0 
showing the most decay and the trapezoidal rule appearing to be the most accurate. 

In figure 2-4, the numerical damping of the various algorithms varies with the product 
ojAt. The idea that we should see less numerical damping at a fixed frequency with a 
decrease in At is illustrated in figure 5-4 which shows the same results comparison as in 
figure 5-3 for a time step of At = 0.001 s. At this time step most algorithms are accurate 
out to the tenth mode. Only backward differences, which due to its first-order accuracy 
is again a poor solution even at very low frequencies, and A“ 2 = 0 (which like backward 
differences is asymptotically annihilating) are worse than this. 

That the other algorithms, with their varying levels of dissipation, have converged to 
the same solution suggests that the remaining error is not due to numerical dissipation. 
Figure 5-5 shows the comparison for four cases with A^ 2 = — | and At = 0.001 s, with a 
varying number of nodes. As the node density is increased, the numerical model is better 
able to resolve the mode shapes associated with the higher frequencies. At n = 800, the 
numerical solution is in agreement with the analytic solution over the full range of the 
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non-dimensional frequency 

Figure 5-3: Power spectra peaks of the response of the free end of the chain for the analytic 
solution and for six variants of the generalized-a method with At = 0.01s, and n = 50. 



non-dimensional frequency 

Figure 5-4: Power spectra peaks of the response of the free end of the chain for the analytic 
solution and for six variants of the generalized-c* method with At = 0.001s, and n = 50. 


81 









Spectral Response Peaks Comparison, X" 2 = -0.5, A t = 0.001 s 



non-dimensional frequency 

Figure 5-5: Power spectra peaks of the response of the free end of the chain for the analytic 
solution and for A ^ 2 = — 5 , At = 0.001 s, with n = 50,200,400,800. Note that the y axis 
scaling has changed from previous power spectra plots. 

analytically computed response. 

These results demonstrate that the ability of the model to accurately resolve high 
frequency response is dependent on temporal and spatial discretizations and on the nu¬ 
merical dissipation for a given algorithm. Given sufficient temporal and spatial resolution, 
all forms of the algorithm appear ultimately capable of accurately calculating the free re¬ 
sponse of the swinging chain. Based on its better accuracy at the larger time step, the 
best choice of algorithm for this problem appears to be the trapezoidal rule. As will 
be demonstrated, however, this may not always be the case, particularly in cases where 
numerical instabilities arise. 

5.1.2 Forced response to an imposed motion 
Two-dimensional simulations 

The second hanging chain problem that we consider is the case studied by Howell [46]. 
In this problem, a 1.75 m long chain is suspended from an actuator which imposes a 
sinusoidally varying horizontal displacement, Q(t), to the top of the chain (see figure 2 - 1 ). 
There is no analytic solution for this problem so we will compare numerical solution 
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Box: n=100, At = 0.01 Trapezoidal: n=100,At = 0.01 Backwards: n=100, At = 0.01 



horiz coordinate (m) 


= -0.7: n=100, At = 0.01 X~ 2 = -0.5: n=100, At = 0.01 = 0.0: n=100, At = 0.01 



horiz coordinate (m) 

Figure 5-6: Snapshots of the chain configuration near the time of expected intersection 
for six variants of the generalized-ct method. 


results to snapshots of the chain configuration derived from experiments conducted by 
Howell. Figure 5-6 shows the configuration of the lower portion of the chain from 3.43 s 
to 3.46 s for six different algorithms, all with n = 100 and At = 0.01 s. Howell observed 
that the free end of the chain intersects the chain above it at approximately 3.4 seconds. 
The box method and trapezoidal rule most closely match this result, with intersection 
occurring by the 3.43 second mark. For the other algorithms, the delay in intersection 
is proportional to the amount of numerical dissipation in the algorithm. The solution 
for backward differences is again the worst; the chain never intersects itself. Likewise 
for == 0, though it comes closer to doing so. For AJ° 2 = —0.7, intersection actually 
happens at 3.47 seconds and for AJ° 2 = —0.5, at 3.5 seconds. 

The situation changes somewhat if we consider the effect of temporal and spatial 
discretization. Figure 5-7 shows the same time points for versions of the box method 
with n = 100,200 and At = 0.01,0.001,0.0001 s. In this case we see that increasing 
the number of nodes does not significantly effect the solution, suggesting that n = 100 is 
adequate to accurately capture the response. An increase in temporal resolution, however, 
from At = 0.01 s to At = 0.001 s, leads to a delay in the crossover to approximately 
3.46 seconds. The result at the even smaller At = 0.0001 s confirms that the solution 
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Box: n=100. At = 0.01 


Box: n=100, At = 0.001 


Box: n=100, At =0.0001 



horiz coordinate (m) 

Box: n=200, At = 0.01 Box: n=200, At = 0.001 Box: n=200, At = 0.0001 



horiz coordinate (m) 


Figure 5-7: Snapshots of the chain configuration near the time of expected intersection 
for the box method with different spatial and temporal discretizations. 

has converged at these smaller time steps. Figure 5-8 shows this same behavior for the 
trapezoidal rule. The only notable difference between trapezoidal rule and box method 
solutions is the better smoothness of the trapezoidal rule solutions at At = 0.01 s. 

Similar results for AJ° 2 = —0.5 are shown in figure 5-9. In this case, the solution at 
At = 0.001 s has not quite converged to the solutions from the trapezoidal rule and the box 
method at the 3.46 second point. The solutions for At = 0.0001 s are in good agreement 
with the similar solutions in figures 5-7 and 5-8. A notable difference in the solutions for 
the various algorithms does appear in the half second (solutions were only run for four 
seconds) following intersection. Both trapezoidal rule and box method solutions required 
significant adaptation of time step to get through the collapse of the lower portion of the 
chain following the crossover. The enhanced stability of solutions with Aj^ 2 = —0.5 allowed 
for a smooth numerical solution in this region, with no or very little adaptation. Without 
experimental verification, however, we cannot say if this numerically more easily obtained 
solution is accurate. 

The basic accuracy of the solutions from all of the algorithms can be verified by com¬ 
parison with Howell’s [46] data for the chain configuration prior to intersection. The data 
points were recovered graphically from digitized versions of the hardcopy plots. Because 
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Trapezoidal: n=100,At = 0.01 Trapezoidal: n=100,At = 0.001 Trapezoidal: n=100,At = 0.0001 


jlj 

0.25 1 ,/\ 

0.25 


°.2 ,i y I 

0.2 

3.43 s 

0.15 

0.15 

-3.44 s 



-3.45 s 

0.1 

0.1 

- 3.46 s 

0.05-- 

0.05 



0 0.1 0.2 0 0.1 0.2 0 0.1 0.2 

horiz coordinate (m) 

Trapezoidal: n=200,At = 0.01 Trapezoidal: n=200, At = 0.001 Trapezoidal: n=200,At = 0.0001 


honz coordinate (m) 


Figure 5-8: Snapshots of the chain configuration near the time of expected intersection 
for the trapezoidal rule with different spatial and temporal discretizations. 


X“ = -0.5: n=100, At = 0.01 X“ = -0.5: n=100, At = 0.001 X“ = -0.5: n=100, At = 0.0001 





0.1 

horiz coordinate (m) 


Figure 5-9: Snapshots of the chain configuration near the time of expected intersection 
for A ^2 = — \ with different spatial and temporal discretizations. 
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t = 2.30 s t = 2.67 s t = 3.07 S 



x/L 

Figure 5-10: Comparison of simulation and experimental results from Howell [46], fig¬ 
ure 5.29. Simulation results are for A^ 2 = —0.5 with At = 0.0001 s and n = 200. 

the original plots did not contain absolute offset information for the points, the experimen¬ 
tal points were aligned with the simulation snapshots by matching the first experimental 
point with the free end of the chain. The comparison for the lower half of the chain is 
shown in figure 5-10. The simulation results are for A“ 2 = —0.5 with At = 0.0001 s and 
n = 200. At this temporal and spatial resolution the solutions from all of the second-order 
accurate algorithms were essentially identical. The results at all three time points show 
good agreement. The comparison at t = 3.07 s improves with a slight adjustment to the 
horizontal offsets that were applied. 

These results are in agreement with the observations drawn from the free response 
problem. At sufficiently small time steps and adequate spatial resolution, all three al¬ 
gorithms: box method, trapezoidal rule, and A^ 2 ~ —0.5, provide accurate solutions. 
Trapezoidal rule is the best choice in terms of the computational costs of accuracy, where 
cost is measured simply in terms of time step. As indicated, however, in regions where the 
solution becomes numerically unstable some numerical dissipation may be necessary to 
obtain a solution. This suggests a trade-off between optimizing the time step for accuracy 
and optimizing the algorithm for stability. 
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- 

10.0 

- 
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Table 5.1: Comparison of the predicted cross-over time and total simulation time before 
failure for three-dimensional simulations of the forced hanging chain. 


Three-dimensional simulations 

In order to further explore these trade-offs, three-dimensional simulations were conducted 
to explore the behavior of the solutions beyond the time when the chain crosses over 
itself. Howell [46] noted that out-of-plane motions of the experimental chain only became 
significant after this point. The simulations were conducted with a small initial out-of- 
plane force applied at the free end to promote the initiation of out-of-plane motion. This 
models the inevitable presence of small disturbing forces which produce instabilities in the 
two-dimensional motion and eventually lead to a fully three-dimensional response. 

Table 5.1 lists the observed time of the chain crossing over itself and the total running 
time (out of a possible ten second simulation) of the simulation before failure. Depending 
on time step, only solutions for —0.4 < A ^ 2 < —0.2 ran for the full ten seconds and resulted 
in an accurate cross-over prediction. At At = 0.01 s, the stable solutions (at Afj = —0.3 
and A ^2 — — 0 . 2 ) were less accurate than the two-dimensional simulations for AJ ° 2 = —0.5 
at this same time step. This is consistent with the observation that as damping increases 
the cross-over time is delayed, until with enough damping it does not occur at all. Also 
consistent with the two-dimensional results is the convergence to an accurate prediction 
of 3.46 s with an increase in temporal resolution to At = 0.001 s. 

The stability of results for A “ 2 = 0.0 and A ^ 2 = 0.1 at At = 0.01 s, but not at 
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time (s) 


Figure 5-11: In-plane and out-of-plane motion of the free end of the hanging chain for the 
stable algorithms with At = 0.01 s. 

At = 0.001 s, illustrates the dependence of the stability on the frequency content of the 
response, the time step, and the damping properties of the algorithm. Because the spectral 
radii in figure 2-4 all initially decrease with the product a;At, at a fixed frequency a decrease 
in At will result in less damping. If the response at that frequency was responsible for the 
instability then the solution at the smaller time step may actually be less stable. 

Figures 5-11 and 5-12 show the in-plane and out-of-plane motion of the free end of 
the chain for all simulations which ran for the full ten seconds. At At = 0.01 s there 
is little consistency between the levels of out-of-plane motion predicted by the different 
algorithms. For the solutions at At = 0.001 s the results for out-of-plane response appear 
roughly equivalent. A trace of the motion of the free end in the horizontal plane for 
At = 0.001 s and AJ° 2 = —0.3 is shown in figure 5-13. The roughly circular whirling 
motion revealed by the trace after the three-dimensional motion is fully developed is the 
type of response that we expect for this problem [68]. 

5.2 Solutions for a full scale mooring 

As a final study of the stability and accuracy of the time integration algorithm we consider 
a full scale mooring with an imposed sinusoidal vertical motion at the top node of the 


88 
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At = 0.001 s: Out-of-Plane Displacement 
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Figure 5-12: In-plane and out-of-plane motion of the free end of the hanging chain for the 
stable algorithms with At = 0.001 s. 



Figure 5-13: Trace of the horizontal motions of the free end of the hanging chain for 
A°° = -0.3 and At = 0.001 s. 












































































mooring. Given the coordinate integration procedure described in appendix E, errors in 
the overall solution will be evident based on the error in the computed coordinates of the 
top node. We can see this if we consider integrals for the top node position in continuous 
form. From equations E.l and E.2 it is clear that we can write those positions as 


x(L,t) = / [1 + e(s, t)] cos <f)(s, t)ds, 

Jo 

(5.9) 

y(L,t) = / [1 + e(s, t)] sin <p(s, t)ds. 

Jo 

(5.10) 


If we had a perfect solution, the dynamic vertical motion at the top, x(L,t ) would al¬ 
ways match the imposed vertical sinusoidal motion and the horizontal motion, y(L,t), 
would always be zero. Ignoring any errors associated with the numerical integration of 
equations 5.9 and 5 . 10 , any deviation away from the ideal solution represents error in the 
computed values of e and (f> along the entire mooring. Thus, comparing the time evolution 
of the computed horizontal displacement of the top node, when the imposed motion is 
purely vertical, provides a simple and convenient estimate of the error associated with a 
particular form of the time integration algorithm. 

The physical characteristics of the trial mooring are the same as for the field exper¬ 
iment, as described in table 4.1 and figure 4-4. For each form of the algorithm under 
consideration, the model was run for 300 seconds of simulation time with a base time 
step of 0.1 seconds. To facilitate comparison with results from a previous version of the 
program, the spatial discretization was uniform over each segment. The flexural stiffness, 
El, of the chain was set to 0.1 N m 2 . The environmental forcing was chosen to simulate 
rather severe conditions with a uniform current of 2.0 m/s and an imposed vertical motion 
with amplitude 2.0 meters and period 8.0 seconds. 

Figure 5-14 shows the computed top node horizontal position for four cases: the orig¬ 
inal box method without any averaging of the coefficient matrices, the box method that 
arises from the generalized-a method with a*, = 0.5, a m = 0.5, 7 = 0.5, the generalized 
trapezoidal rule, and backward differences. For all cases, the values of a*,, a m , and 7 are 
summarized in table 5.2. Note that for both box methods the solution failed before the 
full 300 second run was completed. A solution fails when the nonlinear solver cannot get 
a convergent solution even after adapting At down by a factor of 10 -4 . For the original 
box method without any coefficient averaging this happened at 131 seconds. The box 
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algorithm 



7 

box method w/o coefficient averaging 

- 

- 

- 

box method w/ coefficient averaging 

0.5 

0.5 

0.5 

trapezoidal rule 

0.0 

0.0 

0.5 

backward differences 

0.0 

0.0 

1.0 

A% = 0.0 

0.0 

-0.5 

1.0 

Af 2 = —0.33 (Cornwell and Malkus [20]) 

0.25 

0.0 

0.75 

V 

II 

1 

O 

0.167 

-0.167 

0.833 

0 

II 

■ 82 - 

-0.111 

-0.722 

1.111 


Table 5.2: a*,, a m , and 7 values for the tested algorithms. Solutions for the box method 
without averaging are based on an old version of the program and cannot be obtained 
within the newly developed generalized-a family of algorithms. 

method with coefficient averaging (this is the box method that we can achieve within the 
generalized-a family derived in chapter 2 ) lasted somewhat longer, with failure at 260 
seconds. In addition to those failures, it appears that the solution from the trapezoidal 
rule is beginning to exhibit the same type of behavior starting at around 250 seconds. 
This solution does indeed fail completely at 445 seconds when allowed to proceed beyond 
300 seconds. 

The failure of these three algorithms reinforces two of the important motivations that 
we gave in developing the generalized-a method for the cable dynamics equations. The 
difference in the box method solutions illustrates the improved stability characteristics 
of the algorithm with temporally averaged coefficient matrices. That all three eventually 
failed illustrates the importance of numerical dissipation in improving stability. Figure 5- 
15 illustrates the calculated shear force at the top node for these four trials. With El = 0.1 
for this chain mooring we expect very little shear force. Prior to failure, however, all three 
algorithms developed obvious Crank-Nicolson type noise in the shear force. The solution 
using backward differences has significant error in the calculated horizontal displacement 
and a slightly noticeable drift in the shear force, but remains stable. The numerical 
dissipation associated with backward differences eliminates the Crank-Nicolson noise and 
the solution proceeds with good stability (albeit with relatively poor accuracy). 

Solutions with significantly improved stability and error properties are obtained from 
the one parameter form of the new generalized-a algorithm. Figure 5-16 shows the error 
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Figure 5-14: Calculated horizontal displacement of the top node of the trial mooring for 
the box, trapezoidal, and backward difference algorithms. 


in calculated horizontal displacement for the trial mooring using four different values of 
A^ 2 - Noting the different vertical scales in figures 5-14 and 5-16 it is clear that the drift 
in the calculated horizontal displacement is substantially reduced for all four of these 
cases relative to backward differences. The worst case in figure 5-16 has approximately 
two orders of magnitude less error at 300 seconds. The best case is nearly four orders of 
magnitude better. 

The rate of error growth, defined as the maximum absolute value of the horizontal 
displacement divided by the 300 second simulation time, is plotted (with circles) for a 
number of A“ 2 values in figure 5-17. This error rate is essentially the slope of the trends 
represented by the four curves in figure 5-16. The error is minimized for A^ 2 ~ —0.19. 

Unfortunately, as the additional curves in figure 5-17 illustrate, the optimum value of 
Af/ 2 is highly problem dependent. The first three curves reflect the error growth rate for 
the mooring in the 2.0 m/s current with three different dynamic excitation conditions. The 
second set of three curves shows the error growth rate for the mooring in 0.5 m/s current 
with the same three dynamic excitation conditions. The static shapes of the mooring in 
the two different current conditions are shown in figure 5-18. 




















For the high current configurations, the location of the error minimum does not change 
significantly for excitations with the same period but differing amplitudes. When the 
period of the excitation is changed, the location of the minimum does shift. This behavior 
is consistent with the frequency response of the mooring not changing significantly with 
amplitude of excitation. This contrasts with the low current configurations for which the 
error maxima and minima are shifted most dramatically when the amplitude, not the 
excitation period, changes. 

Such behavior makes it difficult to draw any general conclusions that would aid in 
choosing an appropriate value of A^ 2 for a given problem. We can say that the overall 
level of error appears to be a direct function of the severity of the excitation, as measured 
by the amplitude of the imposed velocity, A for example. The safest choices also seem 
to be AJ° 2 < 0 to avoid the local maxima seen in the low current configuration. 

Additional support for choosing A^ 2 < 0 comes from an examination of the stability 
of the solution as a function of A^ 2 . If we modify our adaptive time-stepping scheme 
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Figure 5-16: Calculated horizontal displacement of the top node of the trial mooring for 
A~ value of 0.0, —0.33, 0.1, and —0.2. 



Figure 5-17: Time rate of growth of the error in the top node horizontal displacement of 
the trial mooring as a function of A^ 2 - 
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Figure 5-18: Static configuration of the trial mooring configurations in 0.5 m/s and 2.0 
m/s uniform current. 

such that it functions like the adaptive relaxation scheme described in appendix C, we 
can determine the largest At that can successfully and consistently be used to propagate 
a solution in time. At each time step, we either increase or decrease At by some small 
factor depending on the success of the solution at that step. Given 

t i+1 = t* + A t\ (5.11) 

if we can successfully solve the nonlinear problem for t t+1 then we increase the time step 

At i+1 = R 2 At\ (5.12) 

and try for the solution at f* +2 = t* +1 + At* +1 . If the solution at f* +1 is unsuccessful, then 
we decrease the time step 

At* 

At* = —, (5.13) 

and try again. R\ and R 2 are constants slightly larger than unity with 1 < R 2 < R\ so 
that a failed step decreases At more than a successful step increases At. For these trials 
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Figure 5-19: Average successfully adapted time step as a function of Af^ over the course 
of a 100 second simulation for the low current configuration with 1.0 m amplitude and 8.0 
second period excitation. 

R 2 = 1.02 and R\ = 1.1. This procedure tends to drive At to an optimum value in a 
relatively small number of time steps. 

Figure 5-19 shows the average (over a 100 second simulation) successfully applied value 
of At as a function of A“ 2 for the low current configuration with 1.0 m amplitude and 8.0 
second period. This configuration was chosen because the simulations with A^ 2 > 0 in the 
latter three curves (the low current configurations) in figure 5-17 required base time steps 
of 0.05 seconds to proceed without constant adaptation 1 . Simulations with A“ 2 < 0 used a 
base time step of 0.1 seconds and proceeded successfully with little or no adaptation. This 
suggested, and figure 5-19 confirms, that the maximum time step value for these cases was 
dependent on Af° 2 . Data from the high current configurations shows a similar trend, with 
the maximum At decreasing sharply for A“ 2 > 0. There is more variability in the data 
for A% < 0, however, as the maximum At is significantly larger than for the low current 
configurations (between 0.5 and 1.0 second) and in each case shows more variability as 
the solution progresses. 

Based on data in figures 5-17 and 5-19 then, we can conclude that a value for A^ 2 


1 Adaptation in those simulations refers to the standard adaptive time-stepping algorithm which reduces 
At by factors of 10 to ensure that the solution remains on the original sample grid. 
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between —0.5 and —0.9 is reasonable in terms of maximizing stability (as measured by the 
largest allowable At) and minimizing the drift error in long time simulations. In this range 
both the error and stability properties appear to be relatively flat and near optimal for most 
of the cases considered. A“ 2 = — 1 should clearly be avoided as it is the box method with 
no dissipation and is prone to the type of failures exemplified in figure 5-15. A trial using a 
value of = —0.98 demonstrated the same failure mechanism after approximately 2250 
seconds of simulation. This suggests that even a small amount of numerical dissipation 
can significantly improve long term stability, but that for guaranteed stability there is a 
nominal level of dissipation which must be provided. A run with AJ° 2 = -0.9 showed no 
signs of Crank-Nicolson noise build-up after 3000 seconds of simulation. 

These results are consistent with the observations gleaned from the hanging chain 
problems, with the additional caveat that in the case of real moorings, solutions with 
A~ -0.5 are significantly more stable than trapezoidal rule solutions. That we might be 
able to use a slightly larger At to achieve the same level of accuracy with the trapezoidal 
rule is no consolation when we cannot in fact get a stable long-time solution at any 
reasonable At. 


5.3 Mesh refinement 

In studying the spatial discretization of a model mooring system there are three important 
factors to consider. At the most basic level we must choose how many nodes to use in 
discretizing each continuous segment of the mooring. The mesh refinement procedure 
described in chapter 3 also requires that we set c w , the weighting factor used in assigning 
the available nodes. Finally, the value of the flexural stiffness, El, for a given segment 
has an important effect on the static solution over that segment. For relatively high El, 
oscillatory solutions for curvature and shear, described in section 3.3, are not typically a 
problem and uniform meshes with relatively low numbers of nodes are generally sufficient. 
For materials with zero El or El just large enough to prevent the singularities associated 
with zero tension, which is the typical situation for chain moorings, these oscillations can 
be quite significant and mesh refinement becomes important. 

For the chain mooring deployed during the field experiment we arbitrarily set the value 
of El to a value of 0.1. Experience has shown that this value is large enough to prevent 
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zero tension singularities in the dynamic solution. In practical terms this is the flexural 
stiffness of a steel wire that is 1.76 mm in diameter. Alternatively, if we take the diameter 
of the chain to be the shaft link diameter and consider that 

d 2 

El ~ EA —, (5.14) 

16 

then our small value of El is 1/6500 smaller than the value of El for a circular rod of 
equivalent axial stiffness. Given that the refined mesh solutions with this value of El are 
satisfactory, it seems reasonable to avoid any question that a larger artificial value of El 
might begin to affect the dynamic solution in a non-negligible way. 

To examine the effect of c w on the static solution on the refined mesh we consider 
two mooring models. The first models the system as a single, continuous shot of chain, 
neglecting the presence of inline instruments. The second models the field experiment 
mooring as it was deployed, with the inline AxPack instruments between shots of chain. 
In both cases the current was uniform over the water column at 0.5 m/s. The static shape 
of the mooring (which is nearly the same for both configurations) for this current profile is 
shown in figure 5-18. For each trial static solution we compare the curvature to a baseline 
solution generated on a uniform mesh with twice as many nodes and El increased to 10.0. 

The static curvature solutions for the continuous chain model are shown in figure 5-20. 
The trial solutions used 162 nodes over the 80.78 m total length of the mooring. The 
solutions on the mesh refined with c w = 10 and Cy, — 50 appear to be a clear improvement 
over the unrefined uniform mesh (cy, = 0) with the same number of nodes and El value. 
To quantify the improvement, the error in curvature for a range of Cy, values is plotted 
in figure 5-21. The error is calculated as the root mean square difference between the 
baseline curvature solution (resampled on the trial solution mesh using cubic splines) and 
the curvature from the trial solution. 

The error is minimized for a value of Cyj & 5. Higher values of c w give too much weight 
to curvature oscillations and produce a mesh which is too coarse in the interior portions of 
the system. This is clearly shown in the top half of the mooring for c w — 50 in figure 5-20, 
where there are now oscillations in the solution where there were none in the uniform mesh 
trial solution. Figure 5-22 shows the mesh density (defined as the number of elements per 
meter) for the same three cases shown in figure 5-20. The solution with the lower weight 
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Figure 5-22: Mesh density after refinement of the all chain mooring. 

is able to distribute sufficient density near the boundary while maintaining a density that 
is not significantly lower than the uniform mesh in the rest of the mooring. The higher 
weight solution devotes many more nodes to the area near the boundaries and as a result 
cannot provide enough density to other areas. 

Figure 5-23 shows the baseline solution, the uniform mesh solution ( c w = 0), and 
two refined mesh solutions for the mooring with inline instruments. The locations of the 
instruments along the mooring are clearly visible as the flat spots in curvature at s ~ 45 m, 
s « 50 m, and s « 57 m. The number of nodes on each of the chain segments was 91 (over 
45.0 m), 18 (over 3.5 m), 36 (over 7.0 m), and 47 (over 23.0 m). Each 0.76 m AxPack was 
modeled using 3 nodes. The baseline solution with uniform mesh had twice the number 
of nodes over each of the chain segments. The relatively larger number of nodes in the 
shorter chain segments reflects the fact that the length of the decay of oscillations in 
the curvature is related more to mesh density than to physical length. This means that 
comparable numbers of nodes must be employed near each segment boundary, regardless 
of the length of the segment. 

For this case there is no striking minimum in the error shown in figure 5-21. A value of 
c w = 20 appears to give the best solution, but values at least out to c w = 100 also appear 
reasonable based on this measure of the error. In looking closely at figure 5-23 for c w = 50, 
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Figure 5-23: Curvature from the static solutions of the mooring with inline instruments. 


however, the sharpness of the plot around the curvature maximum (s « 35) indicates that 
oscillations near the boundaries are being reduced at the expense of an overly coarse mesh 
elsewhere. 


5.4 Comparison with experimental results 

The final phase of the model validation process is a comparison of simulation results to 
data from the full-scale mooring described in chapter 4. For both the two- and three- 
dimensional models we make two types of comparison. In the first we compare time series 
and spectra from individual data sets to verify the ability of the model to accurately 
capture the detailed response of the mooring. In the second comparison we consider 
statistics of the response from all data sets. This analysis provides a check that our 
chosen hydrodynamic coefficients and environmental parameters yield accurate solutions 
over a wide range of forcing conditions. 

The hydrodynamic coefficients for the chain and AxPacks in the validation runs are 
shown in table 5.3. The added mass can be calculated from the added mass coefficients, 
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material 

d (m) 

c dn 

c dt 

Ca n 

a, 

half -inch chain 

0.0495 

0.5 

0.01 

1.0 

0.1 

AxPack 

0.075 

0.8 

0.069 

1.0 

0.5 


Table 5.3: Mass and drag coefficients for the validation simulations. 


C an and C at , according to 


nd 2 

m a„ — Pw ^ D an , 

(5.15) 

_ nd 2 
m at — Pw ^ i 

(5.16) 


where d is the width of a link of chain. Coefficients for the chain are based on experi¬ 
mental results from Gopalkrishnan [37] and previous numerical studies (e.g., [8]). AxPack 
coefficients are approximations based on cylinder and flat plate coefficients. The bottom 
stiffness was set to 100 N/m 2 and the bottom damping ratio to 1.0. The buoy normal drag 
coefficients for the static solutions were 0.5 (in water) and 1.3 (in air). For the purposes 
of the validation, all of these values were chosen because they were physically reasonable 
and produced simulation results that matched experimental results over most data sets. 
Variations on these parameters and schemes for choosing parameters that best match the 
experimental data are studied in detail in chapter 6. 

5.4.1 Two-dimensional model 

For each of the experimental data sets, the effective values for wind and current in the 
two-dimensional plane and the time series of buoy vertical velocity are used as input to the 
model and a time series of mooring response is computed. The procedures for calculating 
these inputs are described in section 4.4. Because of the relatively low currents and 
winds that were observed during the experiment, static solutions for the simulations were 
obtained using the dynamic relaxation procedure described in section C.4. 

Examples of the simulated tension beneath the buoy, along with the corresponding 
experimentally observed values are shown for two cases in figure 5-24. In both cases, the 
agreement between simulation and experiment is excellent. For the 6 December data with 
relatively moderate environmental conditions (approximately 15 knot winds), the mean 
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Figure 5-24: Comparison of experimental and two-dimensional model simulated tension, 
(a) 6 December 1998 at 0800 localtime, (b) 3 January 1999 at 1600 localtime. 

and standard deviation of the simulated tension over the full 200 seconds of simulation time 
(excluding a 10 second initial ramp-up period) were 1503 N and 201 N, respectively. The 
corresponding statistics for the experimental data were 1503 N and 208 N. The statistics 
for the 3 January storm (with near 50 knot winds) also show close agreement: 1611 N 
and 471 N for the simulation compared to 1610 N and 476 N for the experiment. In this 
latter case a few of the tension peaks are higher in the simulation than in the experiment. 
Given the sharpness of these peaks, it is possible that the analog filtering in the buoy 
instrumentation attenuated the experimental signal. 

Figure 5-25(a) shows the tangential acceleration signal recorded by the lowest AxPack 
for the 3 January 1999 storm. For comparison, we calculate the simulated acceleration, 
a(t), at this point based on the tangential velocity, u(t), and the inclination from the 
vertical, 















































01 _I_I_ I _ I _I- 1 - 1 - 1 - 1 - 1 

100 110 120 130 140 150 160 170 180 190 200 

time (s) 


Figure 5-25: Comparison of experimental (a) and two-dimensional model simulated (b) 
acceleration signal at the lowest AxPack for the 3 January 1999 storm event. 

where G is the acceleration due to gravity. This quantity is plotted in figure 5-25(b). After 
time aligning the two signals using the peak observed just before 170 seconds, the results 
look very similar. Based on a comparison of the spectra of the responses (figure 5-26) 
and the experimental standard deviation, 1.54 m/s 2 , and simulation standard deviation, 
1.55 m/s 2 , the level and frequency content of the responses also show excellent agreement. 
The mean of the acceleration (which is an indication of the static tilt of the mooring 
chain at that point) is lower in the experiment than in the simulation (6.54 m/s 2 versus 
8.23 m/s 2 ). This suggests that the AxPack may have been lower along the chain than 
expected. Given the predicted static shape for the mooring under these conditions, any 
error of two to three meters in the position of the lowest AxPack could produce this 
discrepancy in the mean accelerations. 

For a more complete picture of the model performance, we consider the tension statis¬ 
tics for all 119 experimental and simulated data sets. The standard deviation and mean 
of the tension are plotted versus the standard deviation of heave acceleration (a measure 
of the severity of the dynamic forcing) in figure 5-27. The heave statistic is identical for 
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Figure 5-26: Spectral comparison of experimental and two-dimensional model simulated 
acceleration signal at the lowest AxPack for the 3 January 1999 storm event. 


simulation and experiment because the experimental buoy motion is imposed as an input 
for the simulation. Overall, the agreement in the dynamic results (as measured by stan¬ 
dard deviation of tension) is quite good, with nearly exact agreement in low sea states and 
good agreement in higher sea states. The root mean square difference between experiment 
and simulation is 16.1 N. The relative RMS difference, defined as 


e = 


N 


_l 

n 



°~r,exp q 'r,sim 
°Y,exp 


2 

j 


(5.18) 


is 5.8%. 

The simulated mean tensions do not correspond quite as well with experimental re¬ 
sults, but again the trend with sea state appears to be correct. The root mean square 
difference between simulation and experiment, 37.0 N, is less than 16% of the total ob¬ 
served variation in mean tension over the course of the experiment. Given the lack of 
collocation of the wind measurement, the heavy temporal averaging of both wind and 
current, and the assumptions made in projecting wind and current into a two-dimensional 
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Figure 5-27: Comparison of experimental and two-dimensional model simulated tension 
statistics over all 119 data sets, (a) Dynamic response as measured by the standard 
deviation of the tension, (b) Steady-state response as measured by the mean of the 
tension. 


plane, these larger discrepancies in steady-state results are not unexpected. Relative to 
dynamic results for which we have exact knowledge of the forcing (though we are neglect¬ 
ing the horizontal motions of the buoy), we do not have sufficient information to hope for 
an exact comparison. 

Finally, for a frequency domain analog to the time domain comparisons above, we 
consider the errors in the simulated tension spectra. Because the standard deviation is a 
measure of the energy over the entire spectrum, 

S T {u)du, (5.19) 

it is possible for positive and negative errors at different frequency components to effec¬ 
tively cancel in a comparison of standard deviations. A spectral error metric that scales 
similarly to the RMS error in standard deviation, but prevents the cancellation of positive 
and negative errors can be derived by modifying the spectrum from the simulation so that 
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all errors have the same sign, 


‘-’T,sim ( w ) — |<ST,sim( w ) Sj^exp (^01 + >ST,exp( w )- 


(5.20) 


The standard deviation from the discrete form of this modified spectrum with N frequency 
components is 


T,sim 


\ 


ft X/ 

1 


(5.21) 


Analogous to equation 5.18 then, the spectral error over n data sets is simply 


\ 


t n / — _ * 

£ / a T,exp Pr,sim 

n ■ j \ ‘Tr.exp / j 


(5.22) 


Equation 5.22 applied to the full simulation data set produces an error result of 0.068. 
To better understand the magnitude of this error, figures 5-28 and 5-29 show the experi¬ 
mentally observed and simulated tension spectra for the 6 December and 3 January data 
sets. The spectral errors for these two individual cases are 0.040 and 0.074, respectively. 
Visually, the error in these two cases is quite small, indicating that the error value of 0.068 
over the entire data set is quite reasonable. 


5.4.2 Three-dimensional model 

The validation process for the three-dimensional model is similar to that described above 
for the two-dimensional model. Only 60 experimental data sets are available for the 
validation, however, because of the loss of the y accelerometer channel after 27 December. 
For data sets before 27 December we are able to calculate the vertical, horizontal in-plane 
and horizontal out-of-plane velocities of the buoy to use as dynamic inputs into the three- 
dimensional model. Like the two-dimensional simulations, static solutions are obtained 
using the dynamic relaxation procedure. With the current and wind projected into the 
effective plane of the mooring and horizontal motions rotated into in-plane and out-of- 
plane components we can use the same high current static solution as the initial condition 
in all of the dynamic relaxation solutions, regardless of the orientation of the mooring in 
earth reference coordinates. 
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Figure 5-28: Comparison of simulated and experimental tension spectra for the 6 Decem¬ 
ber 1998, 0800 data record. 



Figure 5-29: Comparison of simulated and experimental tension spectra for the 3 Jan¬ 
uary 1999, 1600 data record. 
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Figure 5-30: Comparison of experimental and three-dimensional model simulated tension, 
(a) 6 December 1998 at 0800 localtime, (b) 22 December 1998 at 0800 localtime. 

A comparison of the experimental and three-dimensional model simulated tension be¬ 
neath the buoy is shown in figure 5-30 for the same 6 December data set as in the two- 
dimensional validation and for a storm on 22 December with winds of 35 knots. In both 
cases the results agree well. For the 6 December data set the tension standard devia¬ 
tion from the two-dimensional simulation (201 N) better matches the experimental result 
(208 N). For the 22 December storm the tension standard deviation from a two-dimensional 
simulation is 392 N and the result from the three-dimensional simulation (405 N) is closer 
to the experimental result of 430 N. In both cases, the mean tension is less accurate in 
the three-dimensional simulation than in the corresponding two-dimensional simulation 
(the mean tension in a two-dimensional simulation of the 22 December data is 1571 N). 
Statistics for all of these cases are summarized in table 5.4. 

Tension statistics for all 60 data sets prior to 27 December are plotted in figure 5-31. 
The root mean square difference between experimental and simulated standard deviations 
is 11.2 N. For the mean tensions it is 31.7 N. In the two-dimensional simulations of these 
same 60 data sets, the corresponding differences are 10.7 N and 29.6 N, respectively. 
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data set 

6 December 

22 December 

3 January 

<t t { N) 

T(N) 

MN) 

T(N) 

°r(N) 

T( N) 

experiment 

208 

1503 

430 

1617 

476 

1610 

2D simulation 

202 

1508 

389 

1580 

471 

1615 

3D simulation 

194 

1474 

402 

1552 

- 

- 


Table 5.4: Tension statistics for the comparison data sets. 
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Figure 5-31: Comparison of experimental and three-dimensional model simulated tension 
statistics, (a) standard deviation of tension, (b) mean tension. 

On average then, for the hydrodynamic coefficients and environmental parameters 
chosen for the validation runs, the two-dimensional results are marginally more accu¬ 
rate when compared to experimental data. However, for the purposes of the validation, 
both models appear to accurately simulate the mooring response over a wide range of 
forcing conditions. The primary reason for the different results from the two models is 
that the hydrodynamic coefficients for the simulations were originally chosen to produce 
reasonably accurate results with the two-dimensional numerical model. As discussed in 
section 6.11 the two-dimensional model can often give accurate results with purely vertical 
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input motion, even if the true input motion is three-dimensional, if the drag coefficients 
are adjusted slightly upwards from their true values. In a three-dimensional simulation 
these same values will be slightly too high. This is the situation here where for simplicity 
and consistency we have used the same set of hydrodynamic coefficients for both two- and 
three-dimensional results. 
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Chapter 6 


A Simple Model for Dynamic 
Tension in Catenary Compliant 
Systems 


In this chapter the validated numerical program and data from the field experiment are 
used to develop a simple model to predict dynamic tension in geometrically compliant 
moorings, particularly shallow water oceanographic moorings. Motivated by the strong 
correlation between the tension and acceleration standard deviations in figure 5-27, a 
model is sought that can predict the dynamic tension (as measured by a T , the standard 
deviation of tension) given only very simple inputs. Such a model can ofFer a significant 
reduction in computational cost and provide a framework for the understanding of the 
physics of these systems. While complete time domain simulations have the advantages of 
high accuracy and completeness in terms of resolving the motions and loads throughout 
the mooring, they are computationally expensive. The full set of two-dimensional simula¬ 
tions generated for the program validation in section 5.4.1 took approximately six hours 
to complete on a 533 MHz Alpha LX workstation (119 simulations at approximately three 
minutes per simulation). For analyses requiring long-term statistics of mooring response 
under a wide variety of forcing conditions, as in fatigue studies [39], such an expense can 
be burdensome. In other situations, such as response prediction for offshore floating struc¬ 
tures, a simplified model could eliminate the need for a fully coupled mooring-structure 
interaction model. 
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In the past, analytical formulations for these types of models have been developed for 
the slow drift damping problem. Nakamura et al. [67] used catenary theory to calculate 
the quasi-steady vertical velocity and acceleration along the mooring. By integrating these 
quantities they were able to approximate the dynamic force at the top of the line due to low 
frequency motions in both the horizontal and vertical directions. When investigating the 
role of high frequency dynamics on the damping problem, however, previous investigators 
have relied on numerical simulation [55]. In the development that follows, analytical 
arguments are combined with statistical relationships gleaned from the experimental data 
to develop a model appropriate for wave frequency dynamics. 

6.1 Physical motivation for a simple model 

Previous authors have used a single degree of freedom (SDOF) spring-mass-dashpot system 
to model the dynamic effects in both taut [38] and geometrically compliant catenary 
moorings [34,40]. The equation of motion for the SDOF system shown in figure 6-1 is 

T(t) = Mz(t) + Bz{t) + Kz{t), (6.1) 

where the overdots signify differentiation with respect to time. Reversing the standard 
convention and treating z(t ) as the input and T(t) as the output, the frequency domain 
transfer function, i7(w), for this system is 

H{ u>) = -Mu; 2 + iuB + K. (6.2) 

For a linear time-invariant system, the spectrum of T, St{w), and the spectrum of z, 
S z (lu), are related by 


S t (u>) = \H(u>)\ 2 S z (u). 


(6.3) 
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The spectra of the input velocity, S v (uj), and acceleration S a [w), are related to the input 
displacement spectrum by 


S v (u) = u 2 S z (v), (6.4a) 

S a (u) = u 4 S z (u). (6.4b) 

By substituting equation 6.2 into equation 6.3 and making use of equation 6.4, the spec¬ 
trum for tension can be written as 

S T (u>) = M 2 S a {u) + (B 2 - 2 MK) S v {u) + K 2 S z {lo). (6.5) 


To apply this SDOF spring-mass-dashpot model 
to the data from the SWEX experiment, a nonlin¬ 
ear fitting procedure is used. For each time series 
from the experimental data, spectra of tension 
and heave displacement, velocity, and accelera¬ 
tion are computed. These spectra are then fitted 
to equation 6.5 using a minimization of the spec¬ 
tral error defined by equation 5.22 to determine 



individual values, Mj, Bi, and Ki for that data 
set. On these terms, and elsewhere in this chap¬ 
ter, the subscript i is used to reinforce the idea 


Figure 6-1: An SDOF spring-mass- 
dashpot system. 


that the value in question relates to a single experimental data set. The resulting coeffi¬ 
cients can be plotted against a non-dimensionalized form of the mean tension to observe 
how the coefficients change with the shape of the mooring. The non-dimensionalized mean 
tension, At, is defined as 


T-T 0 
To 


( 6 . 6 ) 


This value serves as a convenient way to represent the amount that the system is pulled 
away from a purely vertical (At = 0) configuration. To is the suspended weight of the 
mooring at slack current: To = wqH, where wq is the wet weight per length of the mooring 
and H is the water depth. 


115 





Figure 6-2: Mass values from each of the 119 spectral fits to equation 6.5. 


Figures 6-2 through 6-4 show the coefficients from the fits to the 119 SWEX data sets. 
The overall quality of each individual fit is quite high. The spectral error over all data 
sets from equation 5.22 is 0.023. The maximum spectral error in any one data set is 0.055 
and 89% of data sets have a spectral error of less than 0.03. There is a significant amount 
of scatter in the aggregated results, however. In spite of the scatter, trends are apparent 
in both the fitted mass and drag coefficients. The mass that participates in the response 
increases with increasing At. This is consistent with additional mooring line being pulled 
off the bottom as At increases. The damping coefficient also increases with mean tension. 
This is a result of both the additional suspended line and the fact that the normal motion 
(and hence normal drag) over the entire mooring increases as the mooring is pulled into a 
more open configuration. There is no apparent trend in the fitted stiffness coefficients. 

The very high scatter in the stiffness is likely due to the difficulty in determining 
a robust value when stiffness effects are relatively small. The scatter in the mass and 
drag coefficients is more interesting, however, as it may well be real. That is, it may 
reflect natural variation in the data that simply looks scattered given the presentation 
as a function of At only. It may also be a reflection of the fact that the model is not 
capturing all of the relevant physics. 
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Figure 6-3: Damping values from each of the 119 spectral fits to equation 6.5. 



Figure 6-4: Stiffness values from each of the 119 spectral fits to equation 6.5. Negative 
values are not physically meaningful; they are an artifact of fitting to data sets with very 
small stiffness effects. 
































Figure 6-5: Static configurations of the simplified SWEX mooring used in the study to 
isolate tension mechanisms. 

From the governing equations (equations A.44 through A.49) the four basic mech¬ 
anisms that produce dynamic tension are inertia, drag, geometric stiffness, and elastic 
stiffness 1 . The spring-mass-dashpot model includes these same mechanisms, but given 
the highly coupled, nonlinear, multiple degree of freedom form of the full model there is 
no particular reason that it should be an accurate SDOF representation of the coupling 
between these mechanisms. To explore these ideas, simulations of a simplified version 
of the SWEX mooring were run with the mooring properties varied so as to isolate the 
various contributions to the dynamic tension. The mooring model consisted of a single 
continuous shot of chain (the AxPacks were removed) in a fixed water depth of 40 m. 

Simulations were run for five levels of non-dimensional mean tension, ranging from 
At = 0.05 to At = 1.0. At each At the static tension at the top of the mooring was 
specified and the static configuration of the mooring was determined using the second of 
the procedures described in section 3.1.1. The static configuration of the mooring at each 
At is shown in figure 6-5. No current was present in the simulations. This procedure was 
used so that At would remain fixed even with variations in the mooring drag coefficients. 
Dynamic excitation was sinusoidal with amplitudes ranging from 0.2 to 2.0 m and periods 


1 Elastic stiffness effects axe negligible in most geometrically compliant systems. 
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variation 

m (kg/m) 

w 0 (N/m) 

c dt 

C dn 

k 

0.01 

31.85 

0.0 

0.0 

mk 

3.73 

31.85 

0.0 

0.0 

tk 

0.01 

31.85 

0.01 

0.0 

nk 

0.01 

31.85 

0.0 

0.5 

mdk 

3.73 

31.85 

0.01 

0.5 

dk 

0.01 

31.85 

0.01 

0.5 


Table 6.1: Variations on the mooring properties used in the simulations to isolate individ¬ 
ual tension mechanisms. Normal and tangential added mass were zero. 


ranging from 4 to 15 seconds. The mooring configurations that were run are shown in 
table 6.1. 

The first four versions can each be used to isolate a single contribution to the dynamic 
tension. For example, with negligible mass, and no drag, the only contribution to the 
dynamic tension in the first variant is stiffness. Because the wet weight cannot be varied 
without changing At, other effects are obtained by subtracting the known stiffness con¬ 
tribution. If Tfc(f) is the dynamic tension record from the simulation with stiffness only 
then the dynamic tension due to mass is 

Tmass = T m k(t) — Tfc(f), (6-7) 

where T m k(t ) is the dynamic tension record from the simulation with both mass and 
stiffness effects present. If cr mass , cr tan , cr nor , and <r s tiff, are the standard deviations of the 
time series of the tension contributions due to mass, tangential drag, normal drag, and 
stiffness, then a convenient way to summarize the effect of each mechanism is to derive 
effective mass, drag, and stiffness coefficients using 


M* = 

°inass 

f 

(6.8) 


O’a 


°tan 

(6.9) 

°dt - 

2 pTrdH (J v j v j 

c* dn = 

^nor 

(6.10) 



1 J TT 3 

^pdHa v \ v \ 

K* = 


(6.11) 


o- z 
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Figure 6-6: (a) Mass, (b) stiffness, (c) tangential drag, and (d) normal drag coefficients 
calculated from simulations with isolated tension contributions. 

Note that for the drag coefficients in particular, these are effective calculated values, 
rather than the actual values assigned to mooring materials for the numerical simulation. 
Standard deviations are used because they are a convenient expression of the amplitude 
of a sinusoidal time series. a a , <J V \ V \, and a z are the standard deviations of the heave 
acceleration, quadratic velocity, and displacement. 

Figure 6-6 shows the four calculated coefficients as a function of At. For each coeffi¬ 
cient type, only simulations that had a symmetric, regular tension response were used to 
calculate coefficients. For the mass coefficients this means that only results for 15 second 
period simulations are used. Simulations with 4 and 8 second periods did not have a 
regular response because of impact loading at the bottom and the lack of damping. For 
the normal drag coefficients only results for amplitudes of 1 m or less and 8 and 15 second 
periods were used. With no inertial forces the tension response at high velocity was not 
symmetric. The full range of simulations were used for the tangential drag and stiffness 
coefficients. 
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All of the coefficients behave roughly as expected. Mass, stiffness, and normal drag 
coefficients all increase roughly linearly with At as additional line is pulled off the bottom. 
The tangential drag coefficient, which at At = 0 is nearly equal to the actual tangential 
drag coefficient used in the simulations, decreases with At. This is because the amount 
of tangential motion along the chain decreases as the chain is pulled into more open 
configurations. 

The coefficients in figure 6-6 represent the behavior of the mooring with little or no 
coupling between the tension mechanisms. The resulting mass and drag coefficients are 
affected by the presence of geometric stiffness, but because stiffness effects are small, 
the results are similar to those that would be obtained if pure isolation were possible. 
Variations mk, mdk, and dk in table 6.1 can be used to calculate mass and drag coefficients 
in the presence of more significant coupling. For these calculations a single effective drag 
coefficient, 


c;, 


^drag 


2 p<£fi r <X„| t ,| 


( 6 . 12 ) 


combining the effects of tangential and normal drag, is used. 

Assuming that the time series of tension for variation mdk (with all effects present) 
can be written as 


T(t) — T mass (t ) + T(jrag(t) + T st ifr(t), (6.13) 

then a drag coefficient in the presence of mass coupling can be calculated by subtracting 
the tension from variation mk (with mass and stiffness) from variation mdk (with mass, 
drag, and stiffness). Likewise, a mass coefficient in the presence of drag coupling can 
be calculated by subtracting variation dk (with drag and stiffness) from variation mdk. 
These results are presented, along with the uncoupled coefficients, in figures 6-7 and 6-8 
for mass and drag, respectively. 

The coupled drag coefficients differ from the more fully isolated results in that they 
represent the drag contribution to tension in the presence of motions which are enhanced 
by mass effects. In the fully isolated case, the drag contribution was calculated in a 
simulation that had no mass. The coupled drag coefficient is calculated by subtracting 
the mass and stiffness contributions to tension from the tension in a simulation with all 
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Figure 6-7: Fully isolated mass coefficient (circles) and mass coefficient when drag is 
present (x). 



At 


Figure 6-8: Fully isolated drag coefficient (circles) and drag coefficient when mass is present 




























effects present. The effects of mass on the motion in this latter simulation are not removed 
and thus the effect of that motion on the drag coefficient is reflected in the final result. 
This same reasoning applies to the coupled mass coefficient as well. 

In both of figures 6-7 and 6-8, the motion effects due to the coupling lead to increases 
in coefficient values. For drag coefficients the presence of mass leads to increased levels of 
motion along the length of the mooring. This increased motion leads in turn to increases 
in the drag forces. Because the calculated drag coefficient is normalized by the motion at 
the top of the mooring only, the increase in the drag contribution to tension is reflected 
by an increase in the drag coefficient. For the coupled mass coefficients, the presence of 
drag restricts the ability of the mooring to deform, in effect increasing the overall stiffness 
of the mooring. To comply with the topside motion then, the amount of mooring line 
pulled off the bottom increases relative to the simulations in which no drag is present. 
This increase in line off the bottom results in a slight increase in the mass. 

The coupling of mass effects into the drag coefficient is clearly the most significant of 
these relationships, particularly at low values of At. This coupling could explain much of 
the scatter that is apparent in the fitted mass and damping coefficients for the experimental 
SWEX data in figures 6-2 and 6-3. Using the the individual coefficients in those figures 
the spring-mass-dashpot model accurately captures the tension response in any single data 
set. This is clear from the good quality of any one of the spectral fits described above. 
However, the coupling between mass and drag means that the coefficients are a function 
both of the steady state configuration and of the excitation frequency and amplitude. 
Thus, when the coefficients are plotted as a function of the configuration (as measured 
by At) they show significant scatter. This scatter, and the underlying dependence on 
both static configuration and input excitation, make it difficult to formulate analytical 
relationships for the coefficients. 

One approach to developing a simple model then is to find representations of the data 
that have low scatter. If the scatter in the data can be minimized, such representations 
could lead to a model that naturally expresses some of the coupling in the system. In such a 
model the coupling is expressed within the form of the model rather than in the individual 
coefficients. This makes mass and drag effects easier to isolate and thus facilitates analytic 
prediction of model mass and drag coefficients. 
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Figure 6-9: Comparison of scatter in the relationship between a T and various function of 
cr a : (a) as a function of cr Q , (b) as a function of Ara a , and (c) as a function of ra a - 

6.2 Development of the simple model 

Figure 6-9 shows three presentations of cr T . In the first, a T is plotted against a a as in 
figure 5-27(a). In the second it is plotted against the product A ra a . The third panel 
presents a T as a function of the product ra a , where r is defined as 


r = 


To' 


(6.14) 


There is a marked reduction in the scatter in this presentation compared to the first panel. 
Motivated by figure 6-9(c) a proposal for the model is 


a T - Mra a + f(a v |„|), 


(6.15) 


where M is a single coefficient that, together with r, expresses the model mass effect for 
any configuration. The simple linear form of the inertia term reflects the trend apparent in 
figure 6-9(c) for values of rc r a < 1.0. This is the inertia dominated regime [40,99]. Beyond 
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Figure 6-10: Comparison of scatter in the relationship between the portion of a T at¬ 
tributable to drag and various functions of u v \ v \. The preliminary value for M, M is 
173.7 kg. (a) As a function of a v \ v \', (b) as a function of Ara v \ v \\ (c) as a function of ra v \ v \. 

this regime drag {f(cr v \ v |)) becomes important and cr T varies away from the straight line 
trend. 

Various forms for f((r v \ v \) can be examined by subtracting a preliminary estimate of 
the inertia contribution from a T . An initial estimate, M', for the value of M is computed 
based on the slope of a line fitted to the data for which ru a < 1.0 in figure 6-9(c). Figure 6- 
10 shows the resulting estimated values for f{o v \ v \) i n the same three presentations as in 
figure 6-9, with cr a replaced by <r v \ v \. 

The scatter in the velocity plots is greater than for the best acceleration case, but the 
relationship in figure 6-10(b), drag as a function of the product Arcr„|„|, appears to have 
the least scatter. It also offers the possibility that a simple linear form can be used to 
model the drag contribution. This form does have the disadvantage that drag disappears 
as At goes to zero. This limitation is addressed more fully in section 6.10. 

With the same type of linear form as the inertia term and the different non-dimensionalized 
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mean tension, the model equation becomes 


a T = M.TO a + -pC d &TdHcr v \ v \. 


(6.16) 


Like M, Cd is a single coefficient for drag in any configuration. The two model coefficients, 
M and Cd can be determined from a linear least squares fit using experimentally observed 
values of a T , a a , er„|„|, r, and At. For n data sets, the formula for the coefficients is 


M 


_\pdHC d _ 



E ( rcr a Arcr„|u|)j 
2=1 

it Tl 2 

E (TO-aArcr^l). E 


E {root 

i =1 


i -1 


i=l 


E 

2=1 

n 

2=1 


(6.17) 


For the 119 data sets from the SWEX experiment, the fitted values are M = 172.8 kg and 
C d = 0.375. 


6.3 Physical interpretation of the simple model 


The variance of tension in the new model is 


o 2 t = (Mr) 2 a 2 + (^pCdArdH^ a 2 ^ + pMTC d ATdH<r a a v \ v \. (6.18) 

Using the linearizing approximation a 2 , , = 3cq) [4] this can be written as 


<4 = (Mt) 2 a 2 + 


3 (^pCdArdH^j cr 2 v + V^pMrC d ATdHa a 




(6.19) 


Neglecting the relatively small covariance between acceleration and velocity to make use 
of the fact that the variance of a sum of independent random variables is the sum of the 
variances, the governing equation for the corresponding physical system is 


T{t ) = Ma(t ) + v 



+ \fZpMCT)dH<J a . 


( 6 . 20 ) 


It is clear from this result that the proposed model can be understood to represent a mass- 
damper system with a linearized damping coefficient that depends on both the quadratic 
drag and inertia. 
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Casting the simple model in the variance form given by equation 6.19 allows for a 
comparison with the terms in the physically motivated SDOF spring-mass-dashpot model. 
For that model, integrating equation 6.5 yields the variance of the tension in data set i as 


= Mfa 2 ai + {Bf - 2 MtKi) < + K\a\. 


( 6 . 21 ) 


Both models represent the dynamic tension as a weighted sum of motion statistics. They 
differ in the coefficient of the velocity term and in the inclusion or absence of the stiffness 
term. 

The qualitative form of the mass term is the same in both models. From figures 6-2 
and 6-7 it is clear that a mass term that grows linearly with non-dimensional mean tension 
is reasonable. Linear fits to either of those results would be of the form Mo + Mi At. From 
a comparison with the model mass term, Mr, it is clear that the implicit assumption in 
the model is that the mass initial value and growth rate are equal. To first order this is a 
reasonable assumption. If the total suspended mass is taken as the mass per length times 
the suspended length, then the r form of the non-dimensionalized mean tension is equal 
to the scope of the mooring, 


T mgL L 
Tq mgH H 


( 6 . 22 ) 


Assuming that the model mass coefficient is equal to the mass per length times the sus¬ 
pended length and that the mass coefficient at Ar = 0, Mo, is equal to mH then for 

-l)=mL (6.23) 


M = M 0 + Mi At 


mH + Mi ( — 
11 


to be true, Mi must equal Mo- 

In the variance form of the model, the coefficient of is 

Z^pC d ArdHj al + V3pMTC d ATdHcr a . (6.24) 

From the time domain form of the simple model, equation 6.20, it can be seen that this 
entire coefficient represents a linearized damping constant. This damping constant can 
be compared to the velocity coefficient Bf — 2MjA,, in equation 6.21 for the spring-mass- 
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Figure 6-11: Total effective damping constant for the experimental spectral data, Bf - 
2 MiKi, and the simple model total damping coefficient from equation 6.24. 

dashpot model. That term represents both a damping and a stiffness effect. 

Figure 6-11 shows the term B'f — 2MiKi for each of the 119 individual fits to the 
SWEX spectral results for the spring-mass-dashpot model along with the total damping 
coefficient for the simple model for each data set calculated from equation 6.24. With 
mass and drag coefficients calculated from a linear fit to the standard deviation form of 
the model, the model total damping coefficient is able to reproduce the nonlinear shape 
and much of the scatter of the spectrally fitted values. With no stiffness, however, the 
simple model does not capture the negative coefficients at low of. For linearized quadratic 
drag in the spring-mass-dashpot model, B 2 oc of [29]. Thus, as the velocity goes to zero 
in this model the intercept of the af coefficient is —2 MK. This term is only important at 
low frequencies and amplitudes where there is little damping. At higher frequencies and 
amplitudes the B 2 term dominates. This higher velocity region is where the simple model 
total damping constant, with its inherent expression of coupling between inertia and drag, 
is accurately reproducing the shape and scatter of the individually fitted values. 

In an undamped spring-mass model, the —2 MK term governs the response near res¬ 
onance. For frequencies above resonance, inertia dominates the response. By neglecting 
this term the simple model is sacrificing accuracy at these lower frequencies. Given the re- 
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versed notions of input and output in the definition of the transfer function in equation 6.2, 
the undamped resonance is defined as the frequency at which infinite wave amplitude pro¬ 
duces zero tension. Thus, neglecting this term is conservative. Additionally, any loss in 
accuracy will be tempered in real situations because there is always some damping present. 
At the very lowest frequencies, the simple model also loses accuracy because it does not 
include a stiffness term like the spring-mass-dashpot model’s K 2 a 2 . This term governs 
the response near zero frequency. 

The relative importance of the two stiffness effects, mass, and damping in the SWEX 
data can be calculated using the coefficients fitted to the spectra of individual data sets in 
section 6.1. The relative magnitude of each of the terms on the right side of equation 6.21 
in comparison with the total tension energy are 


„ Mfal 
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Ji al. ' 

Bfal 


( V B 


f 

J l 


t/ . 


'Ti 


f z . K ? a l 

h < 


(6.25) 

(6.26) 

(6.27) 

(6.28) 


These response fractions are plotted together in figure 6-12. From these fractions it is clear 
that the stiffness term f^ K has little effect over the full range of conditions encountered 
during the SWEX experiment. For low A r the relative magnitude of this term approaches 
20%, but the total dynamic tension in these configurations is relatively low. At higher 
sea states, this term does not contribute significantly to the dynamic tension. /?*, the 
contribution from the stiffness dependent portion of the velocity coefficient is also quite 
small. This explains why the simple model is able to represent the SWEX data without 
any reference to stiffness. The small stiffness effect also explains the high scatter in the 
fitted stiffness coefficients in figure 6-4. 


6.4 Model performance 

To examine the performance of the simple model, three types of analyses are made: 

• Accuracy of the model predicted a T values compared to the experimentally observed 
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Figure 6-12: Portion of the total tension energy attributable to each of the terms in the 
variance form of the spring-mass-dashpot model, equation 6.21. 


values in terms of RMS error, max error, and the number of predictions with error 
less than five percent. 

• Accuracy of tension spectra calculated using a formula derived from the simple 
model. Detailed comparisons are presented for the 3 January 1999 storm data set 
and 6 December 1998 data set. The spectral error over all data sets is also presented. 

• Bootstrap confidence intervals on the fitted model coefficients. 

For the error analysis, the fractional error in the predicted value of a T for data set i is 
defined as 




(6.29) 


(6.30) 


Using these metrics, the RMS error between model fitted and experimentally observed 
values of a T is 2.7%. The maximum error in any one data set is 8.3%. 93% of data sets 
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Figure 6-13: Comparison of model predicted and experimentally observed standard devi¬ 
ation of tension. 

have an error less than 5%. Figure 6-13 shows the model and experimental tension as a 
function of ra a . 

Casting the statistical relationship into the form given by equation 6.19 facilitates the 
prediction of the tension spectrum based on quantities that are easily obtained from an 
input wave spectrum: 

St = (Mr) 2 S a + 3 ^-pCdArdld'j al + V3pMTCd^TdHa a S v . (6.31) 

Comparisons of model predicted spectra calculated using the fitted coefficients and equa¬ 
tion 6.31 with the experimental spectra for the 6 December and 3 January data sets are 
shown in figures 6-14 and 6-15, respectively. For the low sea state case (figure 6-14) the 
response is inertia dominated and the model result agrees well with the experimental spec¬ 
trum across the full range of frequencies. In the high sea state case the basic agreement 
is good, but the model over predicts the spectral peak by 12.5%. Beyond the spectral 
peak, the velocity spectrum falls away quickly while the acceleration spectrum has the 
same basic shape as the tension spectrum; that the predicted tension spectrum is too high 
suggests that the model is over predicting the mass effect for this configuration. 
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Figure 6-14: Comparison of model predicted tension spectra with the experimentally 
observed tension spectrum for the 6 December 1998, 0800 data set. 



Figure 6-15: Comparison of model predicted tension spectra with the experimentally 
observed tension spectrum for the 3 January 1999, 1600 data set. 
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To more fully quantify the spectral performance of the model, the spectral error metric 
defined by equation 5.22 was applied to model predicted spectra for all 119 experimental 
data sets. For each data set, equation 6.31 was used to calculate a model tension spectrum 
for comparison with the experimentally observed tension spectrum. The RMS spectral 
error for all data sets is e* = 0.043. The maximum spectral error in any one data set is 0.10. 
These errors are lower than those for the tension spectra which were calculated from the 
results of the full time domain numerical simulations in the validation in section 5.4.1 (e* 
= 0.068, maximum error of 0.176). They are also not markedly higher than the errors for 
the spectra calculated using the individually fitted coefficients in section 6.1 (e* = 0.023, 
maximum error of 0.055). This is significant because each of those individual fits was 
actually based on a minimization of this same error. Thus, the error for those spectra 
represents a best case which requires the calculation of 119 sets of coefficients. With just 
two parameters for the entire data set, the simple model is able to reproduce the tension 
responses over the entire frequency range with only slightly less accuracy. 

In addition to producing accurate tension results, the fitted coefficients are very robust. 
With their 95% confidence intervals (calculated using the bootstrap method described in 
appendix F) the fitted mass and drag coefficients are 172.8 ± 2.0 kg and 0.375 dt 0.045, 
respectively. The confidence intervals on these fitted coefficients are 1% and 12% of the 
nominal value. These small confidence intervals are important because they support the 
idea that meaningful model coefficients can be calculated using a priori knowledge of 
mooring properties. Large uncertainties on the fitted coefficients would indicate that the 
model was not capturing the scatter in the data. 

These values can be compared to intervals for coefficients for a model based purely 
on the spectrally derived mass and damping coefficients in figures 6-2 and 6-3. In such a 
model, the mass and damping coefficients as functions of At are derived using fits to the 
individual spectrally fitted values. A straightforward example of this type of model is [40] 

M = Mo -I- Mi At, (6.32) 

B = (B 0 + B 1 At)<t v . (6.33) 

Linear fits to these forms yield Mq = 176.5 ± 1.3 kg, Mi = 102.4 ± 21.7 kg, Bo = 
161.3 ± 27.8 kg/s, and B\ = 1255 ± 426 kg/s. At the highest value of At, the confidence 
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intervals for M, and B are 2.8% and 27% of the nominal value. These intervals only take 
into account the uncertainty in the linear fits to the individual coefficients. The actual 
intervals are even larger than this because of the uncertainties in the individual fits that 
are not accounted for in this analysis. These large uncertainties are a result of the highly 
scattered coefficients derived from the individual spectral fits. 

6.5 Model coefficient dependence on physical parameters 

With a validated numerical simulation program it is possible to simulate the entire ex¬ 
perimental data set. This capability permits the calculation of model coefficients for 
parametric variations of the system that was actually deployed. By simulating a large 
number of variations, the dependence of the model coefficients on the system parameters 
can easily be determined. Parameters considered here are the chain normal and tangential 
drag coefficients, Cd n and Cd t , chain normal and tangential added mass coefficients, C a „ 
and C at i anf i bottom stiffness and damping constants. 

The total explored parameter space is shown in table 6.2. In most cases only one 
parameter is varied relative to the baseline case defined in the first line of the table. The 
baseline values are the same as those used for the validation simulations in section 5.4. For 
each set of parameters, the full time domain numerical model is run for the environmental 
conditions in the 119 experimental data sets. Simulations are two-dimensional with only 
vertical (heave) input. Statistics of the tension responses are then computed and a least 
squares fit is used to calculate the model coefficients M and Cd for that parameter set. 
Curves showing the variation in both coefficients while a single parameter is varied are 
shown in figures 6-16 through 6-18. 

Figure 6-16 shows a strong linear dependence of the mass coefficient on both tangential 
and normal added mass. The model drag coefficient also varies with the added mass 
parameters: very slightly with the tangential parameter and a bit more substantially 
for the normal parameter. This dependency indicates that the model form by itself is 
not completely capturing the coupling between inertia and drag. The increase in normal 
motion along the chain that accompanies the increase in mass is causing an increase in 
drag beyond the level accounted for by the coupling term in equation 6.19. Because the 
model input velocity is the same in both cases, the increase in drag must be reflected by 
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variation 

drag 

mass 

bottom 

C'dn 

c dt 

C'a„ 

Cat 

k 

c 

baseline 

0.5 

0.01 

1.0 

0.1 

0.155 

1.0 

1 

0.5 

0.01 

0.0 

0.1 

0.155 

1.0 

2 

0.5 

0.01 

0.2 

0.1 

0.155 

1.0 

3 

0.5 

0.01 

0.5 

0.1 

0.155 

1.0 

4 

0.5 

0.01 

1.5 

0.1 

0.155 

1.0 

5 

0.5 

0.01 

2.0 

0.1 

0.155 

1.0 

6 

0.5 

0.01 

0.0 

0.0 

0.155 

1.0 

7 

0.5 

0.01 

0.5 

0.0 

0.155 

1.0 

8 

0.5 

0.01 

1.0 

0.0 

0.155 

1.0 

9 

0.5 

0.01 

1.0 

0.05 

0.155 

1.0 

10 

0.5 

0.01 

1.0 

0.2 

0.155 

1.0 

11 

0.5 

0.01 

0.0 

0.2 

0.155 

1.0 

12 

0.5 

0.01 

0.0 

0.1 

0.155 

1.0 

13 

0.0 

0.01 

1.0 

0.1 

0.155 

1.0 

14 

0.2 

0.01 

1.0 

0.1 

0.155 

1.0 

15 

0.4 

0.01 

1.0 

0.1 

0.155 

1.0 

16 

0.45 

0.01 

1.0 

0.1 

0.155 

1.0 

17 

0.55 

0.01 

1.0 

0.1 

0.155 

1.0 

18 

0.6 

0.01 

1.0 

0.1 

0.155 

1.0 

19 

0.7 

0.01 

1.0 

0.1 

0.155 

1.0 

20 

0.45 

0.003 

1.0 

0.1 

0.155 

1.0 

21 

0.55 

0.003 

1.0 

0.1 

0.155 

1.0 

22 

0.6 

0.003 

1.0 

0.1 

0.155 

1.0 

23 

0.7 

0.003 

1.0 

0.1 

0.155 

1.0 

24 

0.5 

0.0 

1.0 

0.1 

0.155 

1.0 

25 

0.5 

0.003 

1.0 

0.1 

0.155 

1.0 

26 

0.5 

0.007 

1.0 

0.1 

0.155 

1.0 

27 

0.5 

0.015 

1.0 

0.1 

0.155 

1.0 

28 

0.5 

0.03 

1.0 

0.1 

0.155 

1.0 

29 

0.4 

0.0 

1.0 

0.1 

0.155 

1.0 

29 

0.6 

0.0 

1.0 

0.1 

0.155 

1.0 

31 

0.6 

0.007 

1.0 

0.1 

0.155 

1.0 

32 

0.6 

0.015 

1.0 

0.1 

0.155 

1.0 

33 

0.6 

0.03 

1.0 

0.1 

0.155 

1.0 

34 

0.5 

0.01 

1.0 

0.1 

0.078 

1.0 

35 

0.5 

0.01 

1.0 

0.1 

0.311 

1.0 

36 

0.5 

0.01 

1.0 

0.1 

0.622 

1.0 

37 

0.5 

0.01 

1.0 

0.1 

0.155 

0.0 

38 

0.5 

0.01 

1.0 

0.1 

0.155 

0.5 

39 

0.5 

0.01 

1.0 

0.1 

0.155 

2.0 

40 

0.5 

0.01 

1.0 

0.1 

0.155 

4.0 


Table 6.2: Parameter variations considered in the model coefficient functional dependence 
study. Bold entries indicate a variation in the parameter relative to the baseline value. 





Figure 6-16: Variation of the model mass and drag coefficient with changes to the system 
normal and tangential added mass coefficients. Unless otherwise indicated, all other sys¬ 
tem parameters are at baseline values. For reference, the suspended mass of chain and 
AxPacks at slack current (At = 0) is 161.6 kg. 


an increase in the model drag coefficient. 

A similar effect is evident in figure 6-17 for the model drag coefficients as functions of 
normal and tangential drag parameters. There are clear linear relationships between the 
model drag coefficient and normal and tangential drag. There is also some dependence of 
the mass coefficient on the system drag coefficients. The effect is quite small, however, as 
expected from the earlier analysis of the effect of drag on mass (figure 6-7). 

The dependencies of the model coefficients on the sea bottom parameters are shown 
in figure 6-18. The smallest effect is that of bottom stiffness on model drag coefficient. 
Over a broad range of stiffness levels, the model drag coefficient is nearly constant. There 
is a slight linear increase in drag coefficient with increasing bottom damping. The most 
significant effects of the bottom parameters are on the mass coefficient. As the bottom 
stiffness increases more of the mooring is supported by the bottom, reducing the mass 
of the mooring suspended beneath the buoy. This leads directly to a reduction in the 
model mass coefficient. That the model mass coefficient increases with increasing bottom 
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Figure 6-17: Variation of the model mass and drag coefficient with changes to the system 
normal and tangential drag coefficients. Unless otherwise indicated, all other system 
parameters are at baseline values. For reference, the equivalent tangential drag coefficient 
of suspended chain and AxPacks at At = 0 is 0.015, or nC^ mv = 0.048. 

damping is a result of large accelerations (actually decelerations in this case) of the chain 
near the bottom in the presence of high bottom damping. The resulting increase in inertial 
force is once again reflected in the mass coefficient because of the constant acceleration 
input in the model. 

6.6 Parameter validation using model coefficients 

One potential use of the model coefficients from the parametric studies is to validate 
the choice of system parameters in the time domain simulations. In the validation in 
section 5.4, simulation results were checked against experimental results to ensure that 
the simulation results were correct. Given the numerous parameters in the simulations, 
however, it is conceivable that the right answers could be obtained with several differ¬ 
ent combinations of those parameters. Comparing the model coefficients derived from a 
simulated data set to coefficients derived from the experimental data is one way to check 
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Figure 6-18: Variation of the model mass and drag coefficient with changes to the system 
bottom stiffness and damping parameters. Unless otherwise indicated, all other system 
parameters are at baseline values. 

whether the simulation parameters are actually correct. In a sense this process checks 
not only that the simulation answers are correct, but that they are correct for the right 
reasons. 

Figure 6-19 shows the fitted model mass and drag coefficients for the forty variations 
plus baseline simulation data sets relative to the model coefficients from the experimental 
data set. The small distance between the experimental result and the coefficients for 
variations 17, 18, and 31 suggests that the parameters in those variations more closely 
approximate the true parameters than those in the baseline simulation. Variations 17 
and 18 represent an increase in the normal drag coefficient to 0.55 or 0.6. Variation 31 
increases normal drag to 0.6, but reduces tangential drag to 0.007. Based on this analysis, 
the remaining baseline parameters (for added mass and the sea bottom) all appear to be 
physically reasonable. 

This type of validation cannot be obtained simply by comparing statistical results: 
the RMS difference in tension standard deviation between simulation and experiment was 


138 











































0.9 


0.8 


0.7 


0.6 


•33 


•28 


0.5 


■o 


O 


0.4 


0.3! 

0.2 


•e *36 
•37 


*9*35, 


.. .*19. 

•18 


•27 


.•23 

•3 ’^22 seifs 
•30' 2 i 


25 

^20 

•29 


•11 


0.1 L -13 


0-—‘- 1 - 1 -- 1 -- 1 -1 

155 160 165 170 175 180 185 

mass coefficient 


Figure 6-19: Model mass and drag coefficients for the simulation and experimental data 
sets. Numbers refer to the variation in table 6.2. The experiment result, with confidence 
intervals, is marked by the circle and the dotted box. The baseline simulation result is 
marked by the *. 


5.7% for variations 17 and 18, and 5.8% for the baseline and variation 31. Variation 19, 
the coefficients for which actually fall outside the experimental result confidence intervals, 
also has an RMS error of 5.7%. While this error is minimum for all the simulation data 
sets, analysis of the model coefficients indicates that the normal drag coefficient of 0.7 in 
this variation is too high. 

While this procedure can validate a parameter set as being a reasonable approximation 
to the true parameters, it cannot reveal the true values for those parameters. With an 
exhaustive search of the parameter space, which even for this simple mooring would be 
computationally very expensive 2 , it would be possible to determine the parameter set 
which best matched the experimental results. Given the overlapping confidence intervals 
of all the fitted coefficients, however, the only result that could be accurately reported 
would be the possible ranges of the parameters. 


2 The forty variations in table 6.2, which by no means represent an exhaustive search of the space, took 
approximately eight days to complete on a 533 MHz Alpha LX workstation. 
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6.7 Empirical relationships for the model coefficients 

The strong linear dependence of the model mass and drag coefficients on the system normal 
and tangential added mass and drag parameters suggests the possibility of constructing 
empirical functions which could be used to calculate model coefficients using only the 
known hydrodynamic and material properties of the mooring. Particularly revealing are 
the relationships exemplified in figure 6-17 for the model drag as a function of system 
tangential drag. That the two lines for different normal drag coefficients are separated by 
a constant offset indicates that the model drag coefficient is simply a linear combination 
of the system normal and tangential drag coefficients. 

Ignoring any dependence of the model drag coefficient on system mass or bottom 
parameters, a formula for the model coefficient as a function of system parameters can be 
written as 


C d = {3 dt irC dt + (3 dn C dn . (6.34) 

1 3 dt and (3 dn express the relative weighting of normal and tangential drag in the composite 
model drag coefficient. The factor of n on the tangential term accounts for the definition 
of the tangential drag coefficient based on material circumference, rather than diameter 
as for the normal drag coefficient. The two weights, (3 dt and 0 dn , can be determined by 
a least squares fit to the model results from variations 13 through 33 (in which only the 
drag coefficients were varied) plus the baseline case in table 6.2. The results from the fit 
are 

P dt = 3.79, (6.35a) 

= 0.46. (6.35b) 

Given validated values C dt = 0.01 and C dn = 0.55 from the previous section, these weights 

lead to proportions for tangential and normal effects in the composite model coefficient of 
approximately one-third and two-thirds, respectively. The quality of the fit is quite high. 
The root mean square difference between the actual model drag coefficients and those 
calculated from equation 6.34 with the weights from equation 6.35 is 5.1%. 

The linear dependence of the model mass coefficient on the system added mass co- 
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efficients suggests that a formula similar to equation 6.34 can be derived for mass. The 
dimensional nature of the mass coefficient leads to a more complicated form, however. 
Taking into account the nominal mass that hangs under the buoy in a slack configuration, 
and the rate of increase of suspended chain length with static tension (the slope of the 
line in figure 6-20, </?), a formula for the mass coefficient is 


M — Mtq T 


7T (P 


An, [m + p——C at + Pm n [m + p—C an 


7T d 2 


(6.36) 


The nominal mass, Mt 0 , is defined as the mass plus tangential added mass of all of the 
components hanging beneath the buoy in a slack (purely vertical) configuration. The 
model mass coefficient, M, is a combination of this nominal mass and a weighted sum of 
the virtual tangential and normal mass of additional material that is pulled off the bottom 
as steady state tension increases. The weights, /3 mt and (3 mn , are again determined using 
a least squares fit to simulation results: variations 1 to 12 plus the baseline in this case. 
The fitted weights are 


An, = -0.156, (6.37a) 

An n = 0.102. (6.37b) 


The RMS difference between the actual mass coefficients and the results from equation 6.36 
using these weights is less than one percent. 

To understand the meaning of the various terms in equation 6.36 and the importance 
of the fitted weights, it is useful to consider a uniform mooring in water depth H. Defining 


7T d 2 

T 

7rd 2 


m at — P A Cart 


m a n ~~ P ^ Ca n > 


(6.38) 

(6.39) 


the nominal mass can be written as 


M To = H(m + m at ), 


(6.40) 
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Figure 6-20: Total length of mooring components suspended below the surface buoy as 
a function of static tension. Data are from the two dimensional validation simulations 
presented in section 5.4.1. 

and equation 6.36 becomes 


M = (H + <p(3 mi ){m + m at ) + <pf3 mn (m + m an ). (6.41) 

With this representation, the weights specify or modify the length of mooring material 
that contributes to the tangential or normal mass. This explains why /3 mt is negative. 
The composite model mass is made up of the tangential mass evaluated over some length 
slightly less than the total length of the mooring plus the normal mass evaluated over 
some small length. 


6.8 A priori response prediction 

The primary motivation for developing equations 6.34 and 6.36 is the hope that these 
formulae can be used to calculate the coefficients for a given mooring design based only 
on the known (or estimated) material and hydrodynamic properties of that system. Such 
a facility would permit dynamic response prediction without the costly construction and 
execution of time- or frequency-domain numerical simulations. In cases where the detailed 
information available from such simulations is a necessary part of the design process, 
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response prediction based on the simple model and a priori coefficients could still facilitate 
the early design iteration stages. 


6.8.1 Specifying the steady state tension 

A critical piece of information in the dynamic tension model is the non-dimensional steady 
state tension, r (and At). In early design studies it is probably sufficient to estimate the 
mean tension using catenary formulae. For more refined predictions a static nonlinear 
model, such as the one described in chapter 3 and appendix D, could be run. Easiest of 
all for predicting the response in survivability conditions would be to specify a value for 
r directly. Experience with oceanographic catenary moorings in 40 m or greater water 
depth suggests that a reasonable maximum value for r in similar systems is about 1.3. 

Calculating the mass coefficient (equation 6.36) also requires knowledge of the rate of 
change of the length of the mooring with steady state tension, ip. From the inextensible 
catenary results in appendix G, the rate of increase of suspended length with increasing 
At is 


ip = 


dL 

d.T 


IP 

L ‘ 


(6.42) 


This formula must be employed with some care for non-uniform moorings. More refined 
calculations of <p could be made by running several non-linear static simulations and 
estimating the slope of the resulting (t, L) line, as in figure 6-20. 


6.8.2 Calculating model coefficients 

For the basically uniform all chain experimental mooring, application of equations 6.34 
and 6.36 to calculate model mass and drag coefficients is straightforward - simply input 
the mass and drag properties of the chain. More complicated moorings require some pre¬ 
processing to calculate the input variables for equations 6.34 and 6.36. For a slack mooring 
composed of p segments (chain shots, instrument cages, strongbacks, etc.), equivalent 
normal and tangential drag coefficients are calculated as 




i=\ 


(6.43) 
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di , L{, and C l dn ( are the diameter, length, and drag coefficients of segment i, respectively. 
d is the diameter of the mooring material that includes the grounded portion of the 
mooring. The assumptions behind this approach to averaging the drag coefficients are 
that the mooring is uniform below a certain depth and that the mooring drag coefficients 
can always be characterized by the drag properties of that portion of the mooring that is 
suspended in a slack configuration. The first of these assumptions is not very restrictive. 
Instrumentation is seldom placed below the mud line. The second assumption implies that 
for heavily instrumented large scope moorings at high static tensions, the model prediction 
would be overly conservative. There is no mechanism in the model to account for the fact 
that the long length of ungrounded bottom line in this situation has lower drag coefficients 
than those calculated from the instrumented portion of the mooring. 

For the mass coefficient the process is somewhat easier because there is no averaging 
involved. The nominal mass is calculated from 

p 

M To =^2(m + m at ) i L u (6.44) 

i =1 

where (m + m at ) i is the mass plus tangential added mass per length of segment i. Appro¬ 
priate values for m, C \ n , and C at in equation 6.36 are simply those for the lower uniform 
portion of the mooring. 

6.9 Validation of a priori response prediction 

In order to test the idea that formulae for the model coefficients derived using a data set 
from a single experiment are broadly applicable, three test moorings are considered. The 
first is the shallow water chain catenary mooring from the Coastal Mixing and Optics 
(CMO) experiment, for which experimental results are available. In this case the model 
predictions are compared directly to the experimental results. The remaining test cases 
are contrived examples of an offshore riser in four different configurations: three catenary 
shapes and a lazy wave shape. The lazy wave configuration is the same as that con¬ 
sidered in Larsen’s [61] comparative study of different numerical programs. Because no 
experimental results are available for these cases the model predictions are compared to 
simulation results. 
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6.9.1 CMO mooring 


The central discus mooring of the CMO experiment was an instrumented chain mooring 
deployed in 70 m of water off the northeast coast of the United States from July 1996 
through June 1997 [40]. The central discus buoy contained the same motion package as in 
the SWEX experiment, with a 10,000 pound load cell at the top of the mooring chain. The 
primary difference in the two moorings is the instrument load. In 70 m of water the CMO 
mooring had a nominal virtual mass of approximately 1570 kg of chain and instruments 
suspended below the buoy. The field experiment mooring had approximately 165 kg of 
chain and instruments in 40 m of water. Instrumentation included five vector measuring 
current meters (VMCMs) and four Seacat conductivity and temperature probes. 

The data set from the Coastal Mixing and Optics experiment comprises 634 time series 
of tension and motion. Composite normal and tangential drag coefficients, calculated 
according to equation 6.43, are 


C e d f v = 0.025, (6.45) 

C^ uiv = 0.97. (6.46) 

The nominal mass is 1570 kg and the outside width and mass per length of the bottom 
chain are 0.066 m and 7.98 kg/m, respectively. Applying the weights from equations 6.35 
and 6.37 yields model coefficients of Cd = 0.75 and m = 1553 kg. 

These results are very close to the coefficients calculated from a model fit to the 
564 experimental data sets for which wind data was readily available. From that fit, 
M = 1557 ± 7 kg and Cd = 0.79 ± 0.05. A plot of u T versus ra a for the experimental 

and model results with the a priori calculated coefficients is shown in figure 6-21. Because 

the purpose of the comparison is to validate the a priori coefficients rather than the 
usefulness of the model as an a priori design tool, the model results were calculated using 
the experimental mean tension and a value for <p calculated from 564 static simulations. 
The root mean square difference between the model prediction and experimental result is 
2.1%. This error is the same as that from a comparison of experimental and model results 
with the above mentioned fitted coefficients. 
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Figure 6-21: Comparison of experimental and model predicted a T for the CMO mooring 
using a priori calculated model coefficients. RMS error between experiment and model 
result is 2.1%. 

6.9.2 Catenary riser 

The steel catenary riser (SCR) problem tests the predictive capabilities of the model on 
a problem with a scale typical of offshore energy production systems rather than oceano¬ 
graphic applications. The validation baseline in this case is derived from simulations. The 
simulated system consists of 1500 m of 0.21 m diameter pipe deployed in three different 
configurations. The pipe has a mass per length of 89 kg/m, axial stiffness of 5 x 10 9 N, 
and bending stiffness of 6.6 x 10 3 Nm 2 . Hydrodynamic and bottom parameters in the 
simulation were set to C at = 0, C Qn = 1, Cd t = 0.05, Cd n — 1.0, k = 0.42 and ( = 1.0. 
The simulations were run for vertically imposed motions equal to sea states two through 
nine. Two configurations were run in 600 m water depth with the steady state horizontal 
position of the top node at 1000 m and 1200 m. A third configuration was run in 300 m 
of water with the horizontal position of the top at 1450 m. The current profile in all cases 
was constant at +1.0 m/s (left to right) from the surface to one-third the water depth 
and then decreased linearly to zero at the bottom. The modeled static configurations are 
shown in figure 6-22. These configurations provide working scopes (ratio of suspended 
length to water depth) of approximately 1.1, 1.5, and 3.8. 

For illustrative purposes and because only simulation results are available for compar- 
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Figure 6-22: Static configurations of the catenary riser for the simulation results. 


ison, this example also employs estimates for r and tp based on catenary formulae rather 
than detailed static simulations or experimental results. Thus, the model calculations in 
this case demonstrate the process that a designer might follow in using the simple model 
in the early stages of the design process. Given estimates for the working scopes of 1.1, 
1.5, and 3.8 for the three configurations, equations G.4 and G.6 were used to calculate 
values for r and <p in each configuration. The results of these calculations are shown in 
table 6.3. Because the model mass coefficient, M, depends on p, a different mass coeffi¬ 
cient was calculated for each configuration. Substituting the known pipe properties, the 
weights for the mass coefficient given by equation 6.37, and the different values for p into 
equation 6.36 yields for M of 52527 kg, 52760 kg, and 26574 kg, for the three cases, respec¬ 
tively. Similarly, substituting the pipe drag parameters into equation 6.34 with weights 
given by equation 6.35 yields Cd = 1.055, independent of the configuration. 

Motion statistics for the model were calculated from the same Bretschneider spectrum, 
S(lo), that was used in generating the input time series for the simulation in each sea state, 


roo 

a v\v\ — / u?S{uj)du, 

Jo 


(Tn = 


uj i S{ijS)du 


(6.47) 

(6.48) 
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Table 6.3: 


Configuration scope 


SCR a 1.10 

SCR b 1.50 

SCR c 3.80 

Lazy wave 1.07 


Non-dimensional mean tension and 


r cp (m) 

1.1 545.5 

1.63 400.0 

7.72 78.9 

1.07 290.7 

values for the catenary riser and lazy 


wave riser systems. 


% error in model predicted a T 


sea 

state 

significant 
height (m) 

peak 

period (s) 

config (a) 

config (b) 

config (< 

2 

0.3 

7.5 

-15.1 

-48.6 

-61.5 

3 

0.9 

7.5 

-14.3 

-41.1 

-43.8 

4 

1.9 

8.8 

-11.4 

-29.9 

-27.0 

5 

3.3 

9.7 

-8.3 

-16.3 

-9.0 

6 

5.0 

12.4 

-4.6 

-7.0 

0.9 

7 

7.5 

15.0 

-1.3 

3.7 

12.5 

8 

11.5 

16.4 

3.1 

17.2 

39.2 

8+ 

16.0 

20.0 

0.9 

15.8 

38.7 


Table 6.4: Error in the model predicted a T for the catenary riser using a priori model 
coefficients. Model coefficients were calculated using approximate steady state tension 
results from equations G.4 and G.6. Sea state parameters are based on the North Atlantic 
data from Faltinsen [29], Table 2.3. 


Results of the comparisons are shown in table 6.4 and figure 6-23. The largest relative 
errors occur for the low sea states in all three configurations. At sea states two and three 
in the higher scope configurations the model under predicts the response by 50% or more. 
This represents a much smaller error in the total tension, however, as the static tension is 
very high in these configuration. At higher sea states the agreement between model and 
simulation improves. In the lowest scope case (a) the errors for sea states six and above are 
less than 5%. While the same holds true for the high scope configurations in sea states six 
and seven, the model over predicts <j t by approximately 17% in both sea states eight and 
nine for case (b) and by nearly 40% in case (c). As described in section 6.10, the model 
increasingly over predicts the tension with increasing sea state because the coupling of 
mass into drag actually becomes less at these high sea states. Overall, however, the close 
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sea state 

Figure 6-23: Comparison of simulation and model a T for the catenary riser, (a) Top node 
at x = 1000 m and water depth of 600 m in the simulation; estimated working scope of 
1.1 for model calculations, (b) Top node at x = 1200 m and water depth of 600 m in 
the simulation; estimated working scope of 1.5 for model calculations, (c) Top node at 
x = 1450 m with water depth of 300 m in the simulation; estimated working scope of 3.8 
for the model calculations. 

agreement in both the quantitative and qualitative way in which the model results predict 
the response as a function of sea state suggest that the model can be applied successfully 
to moorings with very different scales than those considered previously. 

6.9.3 Lazy wave riser 

The lazy wave riser problem is based on the configuration described by Larsen [61]. The 
pipe and bottom parameters are the same as for the catenary riser. The water depth is 
355 m and the static position of the top of the riser is x = 350 m, z — 375 m (20 m above 
the water surface). The current flows from right to left (—a; direction) with a constant 
value of 1.0 m/s from the surface to mid-depth and a linear decrease from mid-depth to 
the bottom. The top of the riser is held in place against the current with an applied 
pre-tension. The simulated static configuration is shown in figure 6-24. 
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Figure 6-24: Static configuration of the lazy wave riser for the simulation results. The 
effective working scope for the model predictions was also determined from this result. 

Because of the different shape, the model must be applied with some care in this case. 
The water depth is taken to be the distance from the bottom of the sagged section to the 
top of the riser (above the surface). Likewise, the suspended length is measured from the 
bottom of the sagged section upwards. With these caveats, equations G.4 and G.6 can 
be used to calculate r and y as for the catenary riser. Results of these calculations are 
given in the last line of table 6.3. The calculated model coefficients are M = 36267 kg 
and Cd = 1.055. 

A comparison between simulation and model predicted a T for the same eight sea states 
as for the catenary riser is shown in figure 6-25. The largest relative errors in this case 
occur at the highest sea states, but none of the errors exceed 11%. For sea states six and 
lower the errors are all less than 5%. The good agreement between model and simulation 
in this comparison reinforces the idea that the model is applicable on a range of scales, 
and also suggests that it can be applied to geometrically compliant shapes other than the 
simple catenary mooring. 
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Figure 6-25: Comparison of simulation and model a T for the lazy wave riser. 

6.10 Conditions under which the model breaks down 

While the comparisons above all showed reasonably good agreement between model pre¬ 
dictions and simulation results, there are conditions under which the accuracy of the model 
becomes degraded, such as at the highest sea states in the riser response. To explore these 
conditions the simplified version of the SWEX mooring first introduced in section 6.1 to 
study scatter in the response statistics was subjected to a wide range of forcing conditions. 
In this study, the mooring properties and hydrodynamic parameters were set to their base¬ 
line values. The a priori model coefficients given these properties are M = 156.3 kg and 
Cd = 0.349. Three hundred simulations were run with At values of 0.05, 0.1, 0.2, 0.5, and 
1.0, ten excitation amplitudes ranging from 0.1 m to 2.0 m, and six excitation periods (4, 
6, 8, 10, 12, and 15 seconds). 

Figure 6-26 shows the simulated and model predicted values of a T as a function of a a 
for four values of At. In each case the model prediction agrees reasonably well with the 
simulation for lower values of a a . At all four At values, however, the model over predicts 
the dynamic tension at the highest acceleration levels. The critical acceleration at which 
the model accuracy is significantly degraded increases with At. Thus, both steady state 
configuration and excitation level determine when the model breaks down. In figure 6- 
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Figure 6-26: Simulation and model predicted values for a T in a study using a broad range 
of sinusoidal excitation conditions. 

26(a) for At = 0.05 the model predictions for a a > 2 m/s 2 have relatively large errors. 
For At = 1.0 in figure 6-26(d), only the result at the highest value of a a (approximately 
3.5 m/s 2 ) has a large error. 

The over prediction of the tension is likely due to the presence of the coupling between 
mass and drag in the model. As shown in figure 6-8, the relative importance of the 
coupling on the drag coefficient decreases with increasing At. The model has no way to 
account for this decrease and thus the presence of the coupling leads to an over prediction 
of the tension in severe conditions. That the effect of the coupling should be reduced 
in severe conditions makes sense in that the process whereby an increase in inertially 
induced motions leads to an increase in drag forces should be self-limiting. At some point 
the motion will reach a speed at which quadratic drag, which is proportional to A 2 cj 2 , will 
restrict any additional line motions that might be caused by inertia, which is proportional 
to Alo 2 . The point at which this occurs increases with increasing mean tension because 
as the mooring is pulled open the coupling between inertia and drag is important over a 
broader range of excitation conditions. 
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In addition to this over prediction of tension in severe conditions there are several 
remaining circumstances in which the model cannot accurately predict the dynamic re¬ 
sponse. The two most interesting are both related to elastic stiffness. For moorings with 
inadequate scope, the geometric compliance mechanism can fail and elastic stretching of 
the mooring line becomes important. In these cases, the model would likely under predict 
the tension. The model essentially assumes that there will always be sufficient geometric 
compliance. 

In contrast, for moorings that are basically geometrically compliant, but also relatively 
elastically flexible, the model over predicts the tension. For example, the sea state 8 
simulation result for a T for the higher scope catenary riser problem with the axial stiffness 
reduced by a factor of 100 is a T = 1.05 x 10 5 N. This is nearly three times less than the 
result calculated using the original stiffness. The model has no mechanism to account for 
this reduction. The implicit assumption in ignoring stiffness effects in the development 
of the model is that the mooring line is inextensible. The validity of the inextensibility 
assumption can be checked using the ratio of elastic to catenary stiffness [94], 



where F\ is the horizontal component of tension at the top of the mooring. Results from 
Irvine and Caughey [58] suggest that inextensibility is a reasonable assumption if this 
ratio exceeds 100 - 1000. 

Other failure modes for the model include cases where the mooring is near vertical 
(At « 0) and <j v \ v \ is large or tangential drag effects are substantial. In both of these 
cases, the model mass term will accurately predict the inertial response but the inclusion 
of At in the drag term means that the drag response will be neglected. Also, in cases 
where tangential drag is very large and normal drag is very small, figure 6-6(c) makes clear 
that the total drag coefficient should decrease with At. The model drag coefficient always 
increases. Finally, because geometric stiffness induced dynamic tension is proportional 
to the dynamic length of material off the bottom and steady state forces acting on that 
material, strong bottom currents or very heavy bottom line can increase geometric stiffness 
effects to the point where they become non-negligible. 
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6.11 Effect of horizontal motions on the model coefficients 


All of the simulations used in developing the simple model thus far have used purely ver¬ 
tical input. This approach was based on the derivation of the model using only heave 
statistics. Given the comparisons to simulation results of a wide range of configurations, 
the model is clearly successful at predicting the dynamic tension response to vertical mo¬ 
tions. Furthermore, given the model’s derivation from and success with the experimental 
results from the SWEX and CMO chain catenary moorings, for which the topside motion 
had components in three dimensions, it can be concluded that in these configurations 
the dynamic tension is dominated by the system response to vertical motions. Clearly, 
however, the horizontal motions must produce some contribution to the tension response. 

To explore what effect horizontal motions have on the model coefficients, three-dimensional 
simulations of the experimental mooring were run for the same baseline plus 40 variations 
of the hydrodynamic coefficients listed in table 6.2. Because of the loss of the y accelerom¬ 
eter channel during the experiment there are only 60 available data sets to be simulated 
for each variation. Figure 6-27 shows the fitted model mass and drag coefficients as 
a function of the simulation tangential and normal added mass and drag coefficients for 
both the original two-dimensional vertical only simulations and the new three-dimensional 
simulations. 

The obvious effect of the horizontal motions is an increase in the model drag coefficient 
and a decrease in the model mass coefficient. Based on the convergence to nearly identical 
mass coefficients at zero normal added mass, the change in model mass coefficient appears 
to be due entirely to an effect in the normal direction. In contrast, both tangential and 
normal effects contribute to the increase in model drag coefficient with the addition of 
horizontal motions. Qualitatively, however, there are no significant differences in the coef¬ 
ficients derived from the three-dimensional simulations with both vertical and horizontal 
topside motions. What this indicates is that the model, which uses statistics of vertical 
motion as the only input, can be successfully applied to some systems with horizontal 
motions because the horizontal motions in these systems can be accounted for by changes 
to the model mass and drag coefficients. 

Because of this, figure 6-19, which compared model mass and drag coefficients for the 
vertical motion simulations to coefficients derived from the experimental data set, is not a 
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Figure 6-27: Variation of the model mass and drag coefficient with changes to the system 
normal and tangential added mass and drag coefficients for both vertical and fully three- 
dimensional topside motion input in the simulations. The parameters not under study in 
each panel remain at baseline values. 

good indicator of the true value of the system hydrodynamic coefficients. Rather, it is an 
indicator of the coefficient sets which when used with a vertical motion only simulation 
produce good matches to the experimental results, which are three-dimensional. The best 
choices for drag parameters from that figure were Cd t = 0.01 and C<* n = 0.55. It is now 
clear that these values must be too high. They had to be artificially large to match the 
experimental results to make up for the fact that there was no horizontal motion in the 
simulations. Figure 6-28 shows the mapping of the model mass and drag coefficients from 
the three-dimensional simulations. The best choice in this case is Cd, = 0.003, Cd n = 0.5, 
and added mass coefficients at baseline values. This is variation 25 in table 6.2. 
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Figure 6-28: Model mass and drag coefficients for the three-dimensional simulation and 
experimental data sets. Numbers refer to the variation in table 6.2. The experiment 
result, with confidence intervals, is marked by the circle and the dotted box. The baseline 
simulation result is marked by the *. 

6.12 Horizontal motion effects in very shallow water 

Horizontal motions have a more significant effect on the dynamic tension as the water 
depth decreases. This conclusion became clear during analysis of an experimental data 
set from a 17 m deep National Data Buoy Center (NDBC) test mooring at Duck Pier, 
North Carolina. For that experiment, an instrumented 3-meter discus buoy was deployed 
from July 1997 through January 1998 [88] (due to an instrumentation failure, data is only 
available for the first two months of this period). The buoy contained a six axis motion 
package, current meter, and meteorological sensors. Two load cells and an S4 current 
meter were deployed on the mooring line immediately beneath the buoy. The current 
meter was deployed in the middle of a 7.1 m length consisting of shackle and short shots 
of 1-inch and |-inch chain. The remainder of the mooring line consisted of 41 m of l^-inch 
chain. 

During the analysis of the data from this mooring two things became apparent. First, 
simulation results could not be made to match experimental results without the inclusion 
of the horizontal surge motion at the top of the mooring. There was no choice of hy¬ 
drodynamic parameters which produced an accurate response given only vertical input. 
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Figure 6-29: Experimental and simulated dynamic tension statistics for 126 of the data 
sets from the NDBC Duck mooring. 

Statistical evidence of this inability is shown in figure 6-29. This situation contrasts with 
that for the deeper (40 m) SWEX experiment described above for which the hydrodynamic 
parameters could be increased in conjunction with vertical input to produce a simulation 
that compared well to the three-dimensional experimental result. The second observation 
was that while the dynamic tension model could be fitted to the experimental results the 
fitted drag coefficient was approximately three times greater than that predicted from the 
a priori coefficient prediction procedures outlined in section 6.8 These two observations 
indicate that the presence of horizontal motions in this very shallow water mooring lead 
to a dynamic tension response that is qualitatively different than the response to vertical 
motions that can be characterized by the simple model. 

6.12.1 A model for the dynamic tension response to horizontal motion 

To separate the effects of horizontal and vertical motions as a function of depth, simulations 
of the NDBC mooring were run with horizontal only, vertical only, and combined horizontal 
and vertical input motion at a series of depths from 10 to 40 m. The length of the 
bottom chain was increased at the higher depths so that the touchdown point was always 
away from the anchor. The simulations were two-dimensional, horizontal motion in the 
surge direction only, to minimize the computation time. Tension statistics at six depths 
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are shown in figure 6-30. Three results are presented for each depth: a^ +h , dynamic 
tension in the simulations with both horizontal and vertical input, a”, dynamic tension 
in the simulation with vertical input, and cr£ + er^, the sum of the dynamic tension in 
the simulations with horizontal only and vertical only input. An important observation 
to draw from figure 6-30 is that as the depth increases the results from the vertical only 
and vertical+horizontal simulations appear to converge. This is consistent with the small 
difference between vertical and fully three-dimensional results from the SWEX mooring. A 
second observation is that at lower depths the sum of the dynamic tension from the vertical 
only and horizontal only simulations appear to sum to the results from the simulation with 
both horizontal and vertical input. This is important because it suggests that the effects 
of vertical and horizontal motion on dynamic tension are linearly separable. 

Based on this latter observation then, a modification to the dynamic tension model 
can be proposed as follows: 

a T = Mra a + ~pCdArdHa v ^ + / (horizontal motion statistics, depth). (6.50) 

Figure 6-31(a) shows the dynamic tension statistics as a function of horizontal acceleration 
for the simulations with horizontal only input in 15 m depth. The qualitative similarity 
between this response and the typical response to vertical motions (e.g., figure 5-27(a)) 
suggests a form similar to the model for vertical motions for the horizontal terms, 

<r£ = M h To ax + ^pC%AtS<t V xM . (6.51) 

The superscripts h indicate terms specific to horizontal motion, subscripts x refer to 
statistics of the motion in the horizontal direction, and S is a projected area because it is 
not immediately clear that non-dimensionalizing the drag coefficient using the full water 
depth is appropriate. For that reason it is more convenient to express the model with a 
dimensional drag coefficient as 

<7y = M h Tcr ax + b h At<j V t \ Vt \. (6.52) 

Figure 6-31 (b) which shows a linear trend with quadratic velocity for the initial guess at 
the non-inertial portion of the tension response provides further evidence that this same 
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Figure 6-31: (a) Simulated dynamic tension in the NDBC Duck mooring in 15 m depth 
given horizontal only input motion, (b) Portion of dynamic tension attributable to drag 
with an initial mass estimate based on the slope of the points in (a) with ra ax < 0.8. 

form of model may be appropriate for horizontal motions. 

When the fitted coefficients for this model are mapped over a range of simulation 
parameters, however, the fitted dimensional drag coefficient is insensitive to the value of 
the simulation normal and tangential drag coefficients. For example, at baseline values 
of Ca n = 0.3 and Ca t = 0.003, the fitted dimensional drag coefficient is b h = 1041 kg/m. 
Doubling the normal drag coefficient to 0.6 results in only a slight increase in b h , to 
1080 kg/m. Likewise, increasing the tangential drag coefficient to 0.01 in the simulations 
produces a fitted value for b h of 1044 kg/m. The fitted mass coefficients vary significantly 
with changes to the simulation normal added mass parameter; there does not appear to 
be any sensitivity of the model mass coefficient to tangential added mass. Parameter 
variations were also run with a range of bottom damping and bottom stiffness coefficients. 
The fitted mass and dimensional drag coefficients for some of these variations are listed 
in table 6.5. 

The vertical model form for the horizontal motions, equation 6.52, is only superficially 
appropriate. That the dynamic tension response to purely horizontal motions is not de¬ 
pendent on the drag coefficients or bottom damping suggests that there is no significant 
drag contribution to the tension response. In fact, the only parameter variations that 
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variation 

M h (kg) 

b h (kg/m) 

baseline 

69.3 

1041 

*9 

3 

II 

o 

CTJ 

60.9 

1080 

C dl = 0.01 

69.3 

1044 

C = o.o 

68.0 

1026 

k = 0.056 

67.9 

1180 

= 2.0 

84.9 

1022 

O 

o 

il 

O* 

69.1 

1045 


Table 6.5: Fitted coefficients for the dynamic tension response to horizontal motions using 
the same model form as for vertical motions. Baseline values are C^ n = 0.3, Cd t = 0.003, 
Ca n — 1.0, C at = 0.1, £ = 2.0, k = 0.22. Variations were run for the 15 m depth case. 


produce a significant change in the fitted drag coefficient are changes to the mooring line 
wet weight. Doubling the wet weight (without changing the mass) of all the mooring com¬ 
ponents yields fitted mass and damping coefficients of 97 kg and 2204 kg/m, respectively. 
A better model then is one in which the non-inertial portion of the tension response is 
attributable to a geometric stiffness effect rather than a drag effect. Such a mechanism 
would explain this correlation between the non-inertial portion of the tension and the 
mooring wet weight. A model that makes use of this insight is, 

or£ = M h Tcr ax + k h Ara x . (6.53) 

The form of the stiffness term was chosen because of the strong linearity apparent in 
figure 6-32 for the non-inertial portion of the dynamic tension as a function of A ra x . a x 
is the standard deviation of the surge motion. 

Table 6.6 lists the fitted coefficients for the model described by equation 6.53 for the 
same variations as in table 6.5 plus variations on the wet weight. As expected the stiffness 
coefficient is largely insensitive to changes in any parameter except mooring wet weight. 
There is a slight dependence on the bottom stiffness: a factor of four decrease in k results 
in a twelve percent increase in the fitted stiffness coefficient. The mass coefficient also 
shows a strong dependence on the wet weight. The quality of the fits is not as high and 
the confidence intervals are not as small as for the vertical model applied to the SWEX or 
CMO experimental data. The errors are not unreasonable, however: RMS error between 
the standard deviation of dynamic tension from the baseline simulations and the fitted 
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variation 

M h (kg) 

k h (N/m) 

baseline 

56.5 ± 6.2 

724 ± 100 

C dn = 0.6 

47.0 ± 6.2 

753 ± 100 

t“l 

o 

o 

II 

56.3 dh 6.3 

726 ± 100 

C = o.o 

55.0 ± 6.3 

718 ± 99 

k = 0.056 

55.3 ± 6.0 

815 ± 107 

Can = 2.0 

72.3 ± 6.5 

711 ± 98 

o 

o 

II 

o* 

56.3 ±6.3 

727 ± 100 

wo/wq = 0.5 

51.3 ± 7.7 

493 ± 66 

w 0 /wq = 1.5 

67.6 ± 6.8 

1099 ± 159 

wq/wq = 2.0 

82.4 ± 7.9 

1528 ± 83 


Table 6.6: Fitted coefficients with 95% confidence intervals for the dynamic tension re¬ 
sponse to horizontal motions using the model described by equation 6.53. Baseline values 
for the parameters are given in table 6.5. The wet weight variations are specified as a ratio 
of the specified wet weight to the nominal wet weight in the baseline case. The variation 
is made to the wet weight of all mooring components. 
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Figure 6-33: Simulated and model fitted (equation 6.53) values for the standard deviation 
of tension in response to horizontal input motion. 

results is 25%; 74% of points have an error of less than 10%. Figure 6-33 shows the 
simulated and fitted values for for the baseline case in 15 m depth. 

Equation 6.53 and the strong dependence of the stiffness coefficient on the wet weight 
explain why the horizontal motions in the deeper SWEX mooring contributed so lit¬ 
tle to the total dynamic tension. First, the stiffness term in equation 6.53 scales with 
non-dimensional mean tension, At. The average value of this parameter decreases with 
increasing depth. The maximum value of At during SWEX was 0.17. In the NDBC 
experiment it was 0.99. Secondly, the weight of chain in the SWEX mooring was approx¬ 
imately five times lower than the depth averaged wet weight of components in the NDBC 
mooring. With model stiffness roughly proportional to wet weight, these two factors com¬ 
bine to produce a horizontal motion stiffness effect that is as much as 25 times lower in 
the SWEX mooring than in the NDBC mooring given similar topside motion. 

6.12.2 Parametric dependence of the model coefficients 

A practical benefit to choosing equation 6.53 as the form of the model for horizontal 
motions is that the fitted stiffness coefficients are relatively constant with depth. Figure 6- 
34 shows the fitted stiffness coefficient over a range of wet weight values for depths from 
15 to 40 m. The upward trend in k h is roughly linear with increasing weight, though the 
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Figure 6-34: Fitted stiffness coefficient for the horizontal motion model in 15, 20, 30, and 
40 m water depth. The x-axis is the scaling factor applied to the wet weight of all mooring 
components. 

slope does steepen somewhat as the weight increases. It is difficult to quantify the exact 
relationship between the fitted value of k h and the wet weight in the simulation because of 
the variation in weight over the length of the mooring 3 . Counterintuitively, perhaps, the 
appropriate value is not simply that for the bottom chain. In 15 m depth simulations with 
half the wet weight of the bottom chain, but nominal values for the rest of the mooring, the 
fitted stiffness value was 645 N/m. This value is higher than the 493 N/m stiffness from 
figure 6-34 calculated when all of the mooring component weights were halved, indicating 
that the weights of the components above the bottom chain do effect the stiffness. 

In a simulation with a uniform mooring consisting only of heavy bottom chain with a 
wet weight of 188 N/m, the fitted stiffness coefficient was 850 ±119 N/m. This leads to a 
ratio k h /wq = 4.5. Applying this ratio to the actual mooring with a wet weight averaged 
over a length of 15 m yields a prediction for k h of 693 N/m. This value is slightly low 
compared to the fitted result of 724 N/m in figure 6-34. In 40 m of water, the average wet 
weight increases because of the increased length of heavy bottom chain, and the predicted 
result using the ratio of 4.5 is 788 N/m, which is too high compared to the fitted result 
of 653 ± 173 N/m. All of these values, however, fall within the 95% confidence regions of 

3 The NDBC test mooring consisted of four distinct segments, with heavier segments nearer the bottom. 


164 








Figure 6-35: (a) Dynamic tension response to horizontal motion in the uniform NDBC 
mooring at 15, 25, and 40 m depths, (b) Portion of the dynamic tension attributable to 
stiffness. The initial mass estimate, Mg, is based on a linear fit to the data in (a) for 
Ta ax < 0 . 8 . 

one another. 

To eliminate the difficulty associated with the variation in mooring properties with 
depth and to explore the interdependence of the mass and stiffness coefficients, simulations 
with only horizontal input were run with the uniform version of the NDBC mooring in 
water depths from 15 to 40 m. Because the mass and wet weight properties of a mooring 
line are related through a proportionality constant in most practical situations only the 
mooring mass was varied in these simulations. The wet weight was defined as 

mo = m 5 (l--^^V (6.54) 

Both tangential and normal added mass coefficients were zero. 

With a uniform mooring, the effective wet weight per length and mass per length are 
constant with depth. Under these conditions, there is virtually no depth dependence in 
the fitted coefficients. Figure 6-35 shows the dynamic tension response in the usual way 
for three different depths overlaid upon one another. That the response across depths can 
be plotted meaningfully in the same way as the response at a single depth suggests that a 
fit to the combined data from all depths will yield coefficients that are valid at any depth. 
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Figure 6-36: (a) Mass coefficient fitted to the tension response data in figure 6-35(a) for 
the uniform NDBC mooring at 15, 25, and 40 m depths plus additional results for 20, 30, 
and 35 m depths, (b) Fitted stiffness coefficient for the same data. 

Figure 6-36 shows the fitted mass and stiffness coefficients for the combined data at 15, 
20, 25, 30, 35, and 40 m depths with the mass per length set to 0.5, 1.0, and 2.0 times 
the nominal value of 22 kg/m. The fitted coefficients vary linearly with the mooring line 
mass. The average ratio of model mass coefficient to mooring mass is 2.89 m. The average 
ratio of model stiffness coefficient to mooring wet weight is 4.68. 

The slope of the fitted mass coefficient as a function of mooring mass in figure 6-36 (a) 
has units of length. This length is the amount of mooring chain over which there is an 
inertial response to horizontal motions. Typically, for wave frequency excitation and low 
values of non-dimensional mean tension, only a small region near the touchdown point 
responds with significant acceleration to horizontal motion. The dimensionality of the 
length of this region complicates any attempt to develop a formula for calculating the 
model mass coefficient in an arbitrary system. In the vertical model, the appropriate 
length scale was the water depth. For the horizontal model, the same mass coefficient can 
be applied across depths and thus depth does not provide the appropriate scaling. 

A length scale, l, can be calculated from the ratio of stiffness and inertial effects, 


j _ w 0 a x 
mcr a x 


(6.55) 
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Figure 6-37: Dynamic tension response of the uniform NDBC mooring to purely sinusoidal 
horizontal input motion as a function of depth and excitation period. 

With the wet weight proportional to mg and o ax oc 

l oc (6.56) 

where u p is the peak frequency of the spectrum of horizontal motion. Though this depen¬ 
dence on excitation frequency could cause some of the scatter in the response statistics, 
the variation in excitation frequency over the course of the experiment was not great. 
To verify the dependence on frequency then, simulations were run with purely sinusoidal 
input motion. Figure 6-37 shows the dynamic tension for the uniform mooring in 15 m, 
25 m, and 40 m depth with excitation periods of 5 s, 8 s, and 11s. Excitation amplitude 
ranged from 0.1 m to 1.5 m. For each excitation period, the slope of the response is 
roughly the same, independent of depth. As the excitation period increases, acceleration 
level decreases, and the slope of the response increases. This increasing slope represents 
an increase in the total mass that is needed to keep inertia in balance with stiffness effects. 

The slopes from these results can be compared to values calculated using equation 6.55. 
Given sinusoidal input, cr ax = oj^a x (rather than just being proportional) and the total 
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mass coefficient, ml, becomes 


ml = (6.57) 

This formula yields predicted slopes of 575 kg, 304 kg, and 119 kg for 11 s, 8 s, and 5 s 
excitation periods, respectively. These values compare to average slopes for the results in 
figure 6-37 of 385 kg, 234 kg, and 126 kg. Even within the simulation results in figure 6-37 
the correspondence between frequency and length scale is not exact. The ratio between the 
average slope of the results for 11 s and 5 s cases is 385/126 or 3.06. From equation 6.56 
the expected ratio is (11/5) 2 or 4.84. The conclusion from both comparisons is that 
the calculated length scale becomes too large as the excitation period increases. At the 
highest frequency, where inertia is most dominant, the calculated length scale appears 
to be accurate. As the frequency decreases and stiffness effects begin to dominate, the 
calculated length scale becomes too large. 

With the model mass coefficient written as M h = ml, the length scale l given by equa¬ 
tion 6.55, and the wet weight and mass related according to equation 6.54, the horizontal 
model for a uniform mooring with negligible or no added mass can be written in terms of 
standard deviation of motion only as 

<j£ = WoU x + A Tpk'j , (6.58) 

where (3% is a constant that relates the stiffness coefficient k h to the mooring wet weight. 
Figure 6-38(a) shows the simulation results for 15 m, 25 m, and 40 m depth as in fig¬ 
ure 6-35 along with cr£ calculated from equation 6.58 using a value for /?£ of 4.68 from 
figure 6-36(b). At low values of ra ax the response is inertia dominated and the model 
predictions agree reasonably well with the simulation results. As ra ax increases, however, 
the model prediction becomes larger than the simulation result. This is a consequence of 
overestimating the mass length scale as the response becomes more stiffness dominated. 

While a length scale based on the ratio of stiffness and inertial effects offers a general¬ 
ized procedure for calculating the model mass coefficient, it is not applicable over a broad 
range of conditions. Figure 6-38(b) shows values for a^, calculated from equation 6.53 with 
fitted mass and stiffness coefficients. At high TG ax these results are significantly more ac¬ 
curate than those based on equation 6.58. The process of calculating model coefficients 
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Figure 6-38: (a) Comparison of simulated and model calculated <jy from equation 6.58. 
The simulation results are the same as those shown in figure 6-35. (b) Simulated and model 
calculated a!}, from the original model equation 6.53 with fitted coefficients (M h = 61.4 kg, 
k h = 848 N/m). 

by fitting to experimental or simulation results is much harder to generalize, however, and 
is therefore much less useful in practical applications. 

6.12.3 Practical application of the horizontal model 

Given appropriate coefficients for both the vertical and horizontal models, the separate 
predictions for response to vertical (equation 6.16) and horizontal (equation 6.53 or 6.58) 
motions can be summed to calculate the total dynamic tension response in the presence 
of both vertical and horizontal topside buoy motions. The validity of this approach is 
supported by the linear separability (in a statistical sense) of the response to vertical and 
horizontal motions in figure 6-30. Additional work is required, however, to determine the 
limits of applicability of this approach. 

For the analysis of experimental results, it would be desirable to fit the experimental 
data to a model which combined equations 6.16 and 6.53. The results from such a fit 
would immediately reveal the relative importance of vertical and horizontal effects in a 
given data set. Because of the typically strong correlation between vertical and horizontal 
motion statistics, however, such a fit does not produce reliable results. There are simply 
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too many degrees of freedom in the fit (four) and too little discrimination amongst the 
input parameters. 

6.13 Summary 

While the model for horizontal motions needs to be more fully studied, the overall results 
from the above analyses of simple models for dynamic tension are quite encouraging. For 
the response to pure vertical motion, or for cases with low values of At in which horizontal 
effects can be neglected, the simple model given by equation 6.16 is quite accurate over 
a broad range of conditions. In the analyses above it was applied to a variety of chain 
catenary moorings and steel riser configurations with good success. When combined with 
a validated model for horizontal response effects the range of applicability will be even 
greater. 

As a data analysis tool or as a design tool with a priori predicted coefficients, the 
simplicity of the model is a compelling advantage. In fact, in the latter application, the 
simplicity of the model greatly facilitated the analysis that yielded the rules for a priori 
coefficient prediction. Despite the model’s simplicity, however, it has features which make 
it physically, as well as practically, compelling. In the analysis of these physics, many of 
the important features of the dynamic response of geometrically compliant moorings were 
highlighted: 

• The dependence of the mass term on r and the drag term on At reflects the inertia 
dominated response regime in low to moderate excitation conditions. 

• The presence of the coupling between mass and drag in the model is important in 
the transition between inertia dominated and drag dominated responses. 

• At some excitation level which is dependent on both steady state configuration and 
quadratic velocity, the drag forces overwhelm the inertially induced motions of the 
chain. Under these conditions the coupling term in the model leads to an over 
prediction of the tension. 

• Stiffness effects can typically be neglected at low non-dimensional mean tension, 
except perhaps for very low frequency, large amplitude excitations in which velocity 
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and acceleration are small. Stiffness effects are more important in the response to 
horizontal motions than to vertical motions. 

Unstudied in this chapter is the elastic dominated regime which exists beyond the drag 
regime for cases where the non-dimensional mean tension is high enough to pull all of the 
available line off the bottom. Webster [99] studied this regime in some detail. In these 
cases the system is no longer geometrically compliant and deforms elastically in response 
to dynamic forcing. For the rigid, stiff materials typically used in these systems this can 
be a dangerous regime. This situation can be avoided by designing the mooring with 
sufficient scope given accurate specification of the environmental conditions. 
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Chapter 7 


Bottom Interaction 

In the previous chapter the focus of the analyses was the dynamic tension at the top of the 
mooring. For the most part, the stiffness and damping properties of the bottom played 
little role in determining that response. Previous authors [79] have shown that the bottom 
properties do play a role in governing the response, particularly the bending response, of 
the mooring in the immediate vicinity of the touchdown point (TDP). In this chapter, 
laboratory experiments are used to investigate whether there are excitation conditions 
under which bottom interaction effects do play a role in other aspects of the mooring 
response. Under these conditions, the suitability of the elastic foundation approach in the 
numerical simulations is also investigated. 

7.1 Description of the laboratory experiment 

The laboratory experiments were conducted in the Iselin flume at the Woods Hole Oceano¬ 
graphic Institution. The flume is 20 m long and has a cross-section approximately 1.2 m 
square. It is equipped both with a tow carriage and recirculation pumps, neither of which 
were used for these experiments. The experiments used a section of mooring chain de¬ 
ployed at a fixed position in the flume. Various configurations of the chain were excited 
using a linear servo actuation mechanism. Data was collected from load cells and a digital 
video camera. 
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Figure 7-1: The basic setup for the laboratory experiments. 


7.1.1 Physical layout of the experiment 

The test specimen was a length of j^-inch galvanized steel chain with an outside link 
width of 1.95 cm and a shaft diameter of 0.57 cm. The mass and wet weight of the chain 
were 0.57 kg/m and 4.84 N/m, respectively. The test chain was suspended from the linear 
actuator and run along a bottom platform to an anchor position. The anchor end of the 
chain was held in place using lead weights placed on top of the chain immediately beyond 
the end of the platform. Pretension and excitation levels were constrained so that the 
chain at the anchor end of the platform never lifted off the bottom. Water depth during 
the experiments was 1.1 m. With a bottom platform height of 10 cm, the effective depth 
was 1.0 m. 

A schematic overview of the experiment is shown in figure 7-1. A photograph of the 
physical arrangement of the actuator, lighting, and test specimen is shown in figure 7-2. 
The 10 cm high bottom platform lifts the chain above the tank bottom so that the entire 
chain is in view of the video instrumentation. The platform, a section of wide aluminum 
channel stock, was used with four different surfaces. The basic hard bottom is simply the 
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Figure 7-2: View of the actuator shaft, load cell, test specimen, and lighting arrangement 
looking down the flume from the anchor towards the top of the chain. 
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Figure 7-3: The foam anti-fatigue mat used as a bottom type. 


aluminum covered with black electrical tape to reduce reflectivity. Other bottom types 
were created by placing either a stippled foam or artificial grass mat on top of the tape. 
Photographs of these two surfaces are shown in figures 7-3 and 7-4. The foam material is 
an anti-fatigue standing mat. The artificial grass mat is a green plastic door mat of the 
type commonly used to scrape the bottoms of shoes clean. For a more realistic bottom 
condition, the channel was turned over and filled with sand obtained from West Falmouth 
Harbor. This sand has a relatively uniform grain size of approximately 290 /zm [30]. 

None of these bottoms were soft enough that their stiffness could be easily charac¬ 
terized. The friction coefficient of each bottom was measured by pulling a 90 cm length 
of chain horizontally by hand, at a roughly constant speed, over an approximately 1 m 
length of the bottom. The average of the tension over the duration of the pull was used to 
calculate an estimate of the drag coefficient. Four runs were conducted on each bottom. 
These pull tests were conducted in air; the results are not necessarily directly applicable 
in water, but they do provide a relative comparison of the friction on the different bot¬ 
toms. The results of the four runs, and their average, for each bottom are summarized in 
table 7.1. The hard bottom has approximately one-third less friction than the two mat 
bottoms, which appear to have very similar friction properties within the context of this 
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Figure 7-4: The artificial grass door mat used as a bottom type. 


bottom 

Run 1 

Run 2 

Run 3 

Run 4 

Average 

hard 

0.49 

0.40 

0.45 

0.51 

0.46 

foam 

0.55 

0.78 

0.75 

0.74 

0.71 

grass 

0.61 

0.80 

0.69 

0.75 

0.71 

sand 

1.25 

1.20 

1.39 

1.01 

1.21 


Table 7.1: Friction coefficients, in air, of the various bottom types. 


test. The sand has a high coefficient in this test partly because the chain tends to become 
partially buried over its length as each pull progresses. 

7.1.2 Actuator mechanism 

The actuator is a Parker Hauser HLE-60 with approximately 60 cm of usable throw. The 
actuator is driven through a 4:1 planetary gear box by a Parker Compumotor SM233 
brushless servo motor. The motor is driven by a Parker Compumotor APEX 10 servo 
drive. The test specimen is attached to the actuator carriage via a hardened steel shaft 
that runs through a guide plate at the end of the linear stage. The system is controlled 
by a PC equipped with a Delta Tau PMAC-Lite servo controller card. The PC runs a 
custom designed program which generates the motion profiles, simulating either regular 
or random waves, and downloads them onto the controller card. Once the motion profile 
is started, the process is entirely under the control of the PMAC card which employs a 
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hardware based PID algorithm to command the drive/motor/actuator system. Feedback 
is provided by a 4000 line optical encoder on the motor. Home and limit switches on the 
actuator allow for repeatability to within approximately one millimeter from one run to 
the next. 

7.1.3 Video instrumentation 

One of the significant advantages of working in the laboratory versus working in the field 
is the opportunity to gather data along the whole mooring. The AxPacks on the field 
mooring provide valuable data, but it is impractical to use many more than the three that 
were employed. Also, they only provide relative motion. By using a video system we are 
able to capture the absolute motion of the entire system in a relatively compact and easy 
to interpret data set. 

The video instrumentation system consists of a Pulnix TM-9701 camera and a MuTech 
MV-1500 frame grabber in a 200 MHz Pentium PC equipped with 192 MB of RAM. The 
camera is a progressive scan monochrome CCD camera with electronic shuttering and 
digital 8-bit output via RS-422. It has a resolution of 484 lines and 768 pixels. The camera 
and frame grabber are controlled by a custom written acquisition program that runs on the 
PC. With a relatively simplistic interrupt driven capture algorithm the maximum frame 
rate for full size frames is approximately 15 Hz. With half size frames, which are more 
convenient for processing and storage reasons, the frame rate can be 30 Hz. The frame 
grabber is triggered by a pulse that comes from the servo control computer. 

The post-processing of the imagery is simplified by the use of blacklight and fluorescent 
paint. The mooring chain is painted white and coated with ultraviolet lacquer that fluo¬ 
resces well under black light. During an experimental run all standard lighting is turned 
off and the windows are blacked out. Illumination is provided by six 40 watt blacklight 
fluorescent tubes hanging immediately above the free surface, parallel to the plan view 
of the chain, two 40 watt blacklight fluorescent tubes positioned across the width of the 
tank just above the top of the chain, and a 400 watt theatrical blacklight flood positioned 
behind and above the chain. 
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7.1.4 Force instrumentation 


In addition to the video instrumentation, the model system is instrumented with a small 
waterproof load cell between the actuator arm and the top of the chain. The load cell 
is a Sensotech Model 34 miniature underwater load cell with a 4 - 20 mA output over 
the zero to five pound range of the cell. The current output is dropped across a 5000 
termination resistor to produce a 2 - 10V output signal. This output signal is fed through 
an analog six pole Tschebyscheff anti-aliasing filter with a 20 Hz corner frequency before 
being digitized (100 Hz, 16-bits) and stored. The data capture routine runs on the servo 
control computer as a background process while the motion profile is executing. 

The load cell is attached to the test specimen and the actuator rod using loops of 26 
AWG wire. The top and bottom studs on the load cell have small holes drilled through 
them to accommodate this wiring. The bottom of the actuator shaft also has such a hole. 
A photograph of this arrangement is shown in figure 7-5. The idea behind this attachment 
scheme is to measure the inline tension at the top of the chain. A rigid, vertical connection 
of the load cell to the bottom of the shaft would provide a measurement of the vertical 
component of tension only. 


7.2 Video processing algorithm 


During each experimental run 384 x 242 pixel, 8-bit grayscale video images are captured 
to RAM at 30 Hz. Each image is electronically shuttered at l/60 th of a second. At the 
end of each run, every second frame is written to a compressed disk file, yielding a final 
sample rate of 15 Hz. An example of a single raw image is shown in figure 7-6. Because of 
the fast shuttering the contrast of the image is relatively low. For presentation purposes, 
the image in figure 7-6 was brightened and sharpened using image processing software. 

The raw images are then convolved with a 3 x 3 vertical gradient filter defined as 


1 2 1 

0 0 0 

-1 -2 -1 


(7.1) 


Edges are extracted from the gradient images using a simple threshold. The edge image 
corresponding to the raw image in figure 7-6 is shown in figure 7-7. At any given horizontal 
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Figure 7-5: View of the actuator shaft, load cell, test specimen, and bottom platform 
through the glass wall of the flume. 



Figure 7-6: Example of a raw image from the video capture system. 
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Figure 7-7: Edges extracted from the raw image in figure 7-6. 


position, the vertical centerline of the chain at that position is calculated as the median 
location of all points along a vertical line. This procedure reduces the edge image to an 
image with no more than one pixel illuminated per horizontal coordinate. These pixels 
are turned into a line through a simple connection of adjacent points. The result of this 
final processing stage on the example image is shown in figure 7-8. 

7.3 Mooring dynamics in the touchdown region 

The initial series of experiments were all conducted on the basic hard bottom described 
in section 7.1.1. Each experimental run lasted twenty seconds, with a two second linear 
ramp of the excitation amplitude at the beginning and end. Excitation amplitudes were 
0.1, 0.15, 0.2, or 0.25 m. Excitation periods were 1.25, 1.5, 2.0, and 3.0 seconds. These 16 
excitation conditions were run at non-dimensional mean tensions, At, of approximately 
0.16, 0.37, and 0.80. At is defined by equation 6.6, 

For reasons of brevity, only results for 0.25 m excitation amplitude and the highest 
and lowest mean tensions are presented. All of the different qualitative dynamic features 
are evident in this subset of the results. Time series of tension for the highest At and 
lowest At values are shown in figures 7-9 and 7-10, respectively, for excitation periods 3.0, 
2.0, and 1.25 seconds. In both cases there is a marked difference in the tension response 
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Figure 7-8: Line representing the center of the model chain extracted from the edge image 
in figure 7-7. 



Figure 7-9: Tension time series for the hard bottom at At ~ 0.80 for excitation amplitude 
0.25 m and excitation periods (a) 3.0 s, (b) 2.0 s, and (c) 1.25 s. 
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Figure 7-10: Tension time series for the hard bottom at At & 0.16 for excitation amplitude 
0.25 m and excitation periods (a) 3.0 s, (b) 2.0 s, and (c) 1.25 s. 

between the slowest and fastest excitation levels. For the 3 second excitation cases, the 
response is roughly sinusoidal, matching the regular input motion. As the excitation 
period decreases, however, the tension response becomes more and more asymmetrical. 

To more fully understand what is happening in the high frequency excitation cases it is 
instructive to consider the motion and tension of the chain over a single cycle. Figure 7-11 
shows the positions of the chain extracted from the video and the corresponding tension 
record for a single cycle of motion starting at 13 seconds for At « 0.80. The top left panel 
shows the chain positions while the motion of the top of the chain is upwards (vertical 
velocity greater than zero) with the starting position drawn in bold. The top right panel 
shows the chain positions during the downward motion, with the first downward position 
drawn in bold. In the tension plot, time points marked with circles correspond to the 
timing of the upward moving position snapshots; squares correspond to downward moving 
snapshots. Starting from the lowest point in the motion, the tension very gradually 
increases until approximately 13.2 seconds at which point it increases very rapidly. The 
tension remains relatively high for approximately 0.15 seconds before falling gradually 
until 13.6 seconds. After that point the tension increases very slowly for the remaining 
0.6 seconds (nearly half) of the cycle. 

At the beginning of the cycle the input velocity is zero and the chain top is at its 
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Figure 7-11: Chain response on the hard bottom over one cycle at 1.25 s excitation period, 
0.25 m excitation amplitude, and At « 0.80 The bold line in each of the top panels marks 
the first profile of that panel. The arrow indicates the direction of motion of the top of 
the chain. In the tension plot, circles correspond to the time points of the upward moving 
profiles, squares to downward moving profiles. The dashed line is the static tension level. 
Dotted vertical lines mark the T p / 4, T p / 2, and 3T p /4 points. 

lowest point. As the chain moves upwards, drag increases as velocity increases. The large 
jump in tension just after 13.2 seconds is due not to drag, however, but to a snap load 
that occurs when the slack, grounded chain suddenly retensions. This phenomenon can 
be seen clearly in close-up video of the touchdown region in figure 7-12. This imagery 
was actually taken for a slightly different case (artificial grass bottom which was held in 
place by a light coating of sand), but the features and timing are nearly the same as in 
the hard bottom case. As the chain moves downward in the moments preceding the cycle 
under consideration, the chain that is being grounded is slack. By the 13.14 second image 
the input motion has started moving upwards again and it is clear that the slack in the 
grounded chain is beginning to be pulled out. When it is fully pulled out, the tension spike 
occurs. Drag keeps the tension relatively high for a time because the bulk of the chain 
is moving very fast, as evidenced by the large separation between profiles in the upward 
moving panel in figure 7-11. As the chain slows in its upward motion, drag decreases and 
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Figure 7-12: Closeup view of the touchdown region showing a sequence in which the chain 
is laid down with slack and then pulled taut. For practical reasons, the bottom in this 
case was the artificial grass mat with a light coating of sand to hold it in place. As will be 
shown in section 7.4, the results for this bottom are nearly identical to those on the hard 
bottom. 











tension decreases. 

At the transition from upwards to downwards motion near 13.6 seconds, the velocity is 
zero and the displacement and acceleration have maximum magnitude with opposite signs. 
Given an acceleration of nearly two-thirds the acceleration due to gravity, the inertial effect 
greatly reduces the increased tension attributable to the weight of the additional line that 
is pulled off the bottom. Thus, with little drag, the tension at this mid-point in the cycle 
is very low. The tension remains low after this point because at this point the chain that 
is being grounded is laid down slack. With no tension at the bottom of the chain, the 
curvature near the top of the chain reverses as the downward motion progresses. With 
the chain top more horizontal and the mid-section of the chain moving relatively slowly 
(as evidenced by the close spacing of profiles in the downward profiles of figure 7-11) due 
to this curvature reversal, there is little dynamic contribution to the tension during this 
part of the cycle. 

Both tension discontinuities, the spike just after 13.2 seconds, and the slack at the 
touchdown point starting at 13.6 seconds, are the result of a shock in the tension. Using 
an analytical result for the interaction of string and bridge in a sitar by Burridge et al. [12], 
Triantayllou et al. [94] predicted that for the cable bottom interaction problem, shocks 
will occur when the velocity of the TDP exceeds the speed of transverse waves in the 
cable. Essentially, the transverse wave speed governs the ability of the mooring line to 
comply geometrically with a smooth rolling and unrolling motion. When the touchdown 
point moves faster than this speed during loading (upward motion) snap loads occur. A 
shock during unloading (downward motion) produces a slack condition at the touchdown 
point. Both of these conditions can be seen quantitatively in the experimental results. 

Following a result from Burridge and Keller [11], 
this shock criterion can be derived by considering the 
integral form of the momentum equation for the sit¬ 
uation diagrammed in figure 7-13, 

j t J b m—ds = T(b , t) - T(a, t) + F + f f(s)ds. 

(7.2) 

The limits a and b define a small region that contains 



Figure 7-13: Definitions for the 
derivation of the shock criterion. 
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the instantaneous TDP which is located at s = so(t). T is the line tension, x(s,t) = 
[:c(s,i),z(s,t)] T describes the position of a point on the line, f is the force density due 
to weight and buoyancy, and F is a reaction force exerted on the line by the bottom at 
the TDP. Using Leibnitz’s rule [43] the integral on the left side of equation 7.2 can be 
evaluated as 
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(7.3) 


Assuming that the line comes instantaneously to rest after being grounded, the velocity 
of the line immediately to the left of the TDP is zero and 



(s 0 ,t) = 0. 


(7.4) 


Furthermore, letting a and b approach so from below and above, integral terms go to zero 
and equation 7.2 becomes 

-m^-^x(s 0 ,t) = T (s^)-T (s^) + F. (7.5) 

Finally, noting that there is no vertical component of tension to the left of the TDP, and 
that for small vertical displacements 


cos (f) = 


x PS s, 
dz 


ds ’ 


(7.6) 

(7.7) 


the force balance in the vertical direction is 


dx o d . . m d . . „ 

■ m rfrs 2(xo ’ !)=r °a; z(l “’' ) + F ' 


(7.8) 


Because z(xo(t),t) is zero for all t, the total derivative of z at the TDP must also be 
zero, 


d , . dx o d , . 


(7.9) 
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Substituting this geometric constraint into equation 7.8 yields 


m 





(7.10) 


From geometric considerations and the assumption that F is an upwards directed reaction 
force, all of the terms in this equation are positive or zero. This leads to two possible 
scenarios: 


m 


m 



>T 0 =» 

d 

—z(x 0 ,t) > 0, 

F>0, 

(7.11) 

<T 0 => 

0 

—z( x o,t) = 0, 

F = 0. 

(7.12) 


In the second scenario, the line leaves the bottom tangentially and there is no impact 
force. The first scenario is the case in which a tension discontinuity forms. The condition 
for this case can be re-written as 



(7.13) 


The quantity on the right in equation 7.13 is the transverse wave speed in the line. When 
this condition is true, there is an impact from the bottom and the line does not leave 
the ground tangentially. This impact force and the loss of tangency introduce the tension 
discontinuity that is the most obvious consequence of the shock. It is important to note 
that while the impact force in this derivation is not itself evident in the topside tension 
record, it does have direct implications for the numerical simulations, as discussed in 
section 7.5. 

Both the transverse wave speed and the TDP speed can be calculated for the experi¬ 
mental results. For the wave speed, To can be estimated by the horizontal component of 
the top tension. The TDP speed is calculated by numerically differentiating the horizontal 
TDP coordinates, xq, extracted from consecutive chain profiles. Figure 7-14 shows these 
two results for the same high frequency, high At case as in figure 7-11. The exceedance of 
the shock criterion, equation 7.13, is clear at both the 13.2 and 13.6 second time points. 

The utility of this criterion in predicting these tension discontinuities is further sup¬ 
ported by the data from the the 2.0 and 3.0 second excitation period cases. In figure 7-15 
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Figure 7-14: Transverse wave (solid line) and TDP (dashed line) speed over one cycle at 
1.25 s excitation period, 0.25 m excitation amplitude, and Ar « 0.80. Circles and squares 
indicate upwards and downwards input motion as in figure 7-11. 

for T p = 2.0 s there is no snap load during the upwards motion but the tension does exhibit 
the slacking response during a portion of the unloading half of the cycle. Correspondingly, 
in figure 7-16 the TDP speed exceeds the estimated wave speed during unloading, but not 
during loading. Note that with slack in the grounded chain, the horizontal component of 
the top tension overestimates To and the TDP speed likely exceeds the wave speed for 
some length of time beyond the brief exceedance shown in figure 7-16. This estimate is 
valid up to the point of the criterion being met, making it useful for the predictive purpose 
shown here, but is not accurate once the tension discontinuity has formed. The response 
in this case also differs from the T p = 1.25 s case because the lower frequency excitation 
leads to a basic tension response that is not simply drag dominated, with weight and 
inertia effects largely canceling one another. 

For 3.0 s period excitation, neither snapping nor slacking behavior is evident in figure 7- 
17. This is expected as the TDP speed in figure 7-18 never exceeds the transverse wave 
speed. At this lowest frequency the tension response is dominated by geometric stiffness. 
The phase of the tension in figure 7-17 very nearly matches the phase of the displacement 
of the chain top. 

The results for 1.25 s excitation period and Ar ss 0.16 (the lowest non-dimensional 
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Figure 7-15: Chain response on the hard bottom over one cycle at 2.0 s excitation period, 
0.25 m excitation amplitude, and At « 0.80. 



Figure 7-16: Transverse wave and TDP speed over one cycle at 2.0 s excitation period, 
0.25 m excitation amplitude, and At « 0.80. 
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Figure 7-17: Chain response on the hard bottom over one cycle at 3.0 s excitation period, 
0.25 m excitation amplitude, and At 0.80. 












Figure 7-19: Chain response on the hard bottom over one cycle at 1.25 s excitation period, 
0.25 m excitation amplitude, and At « 0.16. 

mean tension) are shown in figures 7-19 and 7-20. Qualitatively, the response in figure 7- 
19 is similar to that for the At « 0.80 case in figure 7-11. The onset of the snap load 
is delayed relative to that case because the higher initial curvature of the low tension 
configuration at its lowest point results in the TDP speed reaching its maximum more 
slowly. The slack discontinuity occurs at the same time in the two cases because that 
shock is more dependent on a low wave speed than on a high TDP speed and the phase 
of the wave speed is similar in the two cases. 

7.4 Effect of bottom conditions on mooring response 

7.4.1 Artificial bottoms 

In addition to the hard bottom tests described above, tests were run on the artificial 
bottom types described in section 7.1.1. The artificial mats have higher friction than the 
hard bottom and some unquantified differences in their stiffness properties. Based on 
the results, however, these properties are only weakly relevant (a conclusion supported 
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Figure 7-20: Transverse wave and TDP speed over one cycle at 1.25 s excitation period, 
0.25 m excitation amplitude, and At « 0.16. 

by the full scale experimental and simulation results in chapter 6) or the differences were 
not significant enough to produce a marked change in the response. Figures 7-21 and 
7-22 show the tension time series for the runs on foam and artificial grass bottoms at 
At cs 0.80 with excitation period 1.25 seconds and excitation amplitude 25 cm. There 
are no significant differences between these results and those shown in figure 7-9 for the 
hard bottom case. The mean values are slightly different due to the added height of the 
bottom mats and the accompanying small variations in the shape of the chain. 

7.4.2 Sand bottom 

A more interesting response was observed in the runs on the sand bottom. Like the 
artificial bottoms, the tension records for these runs do not look markedly different from 
those obtained on the hard bottom. The interesting feature of the response on sand is 
the trenching and digging action of the cycling chain. Because of this action, the chain 
was often below the plane of the bottom and thus was not visible to the camera. For this 
reason, the standard high speed video and associated processing were not used for runs on 
sand. Instead, the camera was repositioned to look down at an angle on the touchdown 
region (this is the position from which the closeup video in figure 7-12 was taken). To 
document the trenching behavior, time lapse video was then taken every ten seconds over 
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Figure 7-22: Tension time series for the grass bottom at At ~ 0.80 for excitation amplitude 
0.25 m and excitation periods (a) 3.0 s, (b) 2.0 s, and (c) 1.25 s. 
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Figure 7-23: Tension time series of the initial twenty seconds for the sand bottom with 
At « 0.16, excitation amplitude 0.25 m, and excitation periods (a) 3.0 s, (b) 2.0 s, and 
(c) 1.25 s. 


the course of two consecutive three minute runs. The sand was restored to its original, 
flat condition after every six minute run (between each change in excitation period or 
non-dimensional mean tension). Because this process was more time consuming than the 
runs on the artificial bottoms, only 25 cm excitation amplitude cases were performed. 
Tension data was captured as before at 100 Hz. 

Figure 7-23 shows the first twenty seconds of the tension record for the runs on sand at 
At ~ 0.16. In the high At runs, the lowering of the bottom as the trench deepened over 
time, and the subsequent rise in steady state tension, led to tension spikes in the 1.25 s 
excitation period case which were clipped in the data acquisition system (over 5 lbs). For 
this reason, the low At runs are used to facilitate a direct comparison with the results 
already presented for the hard bottom runs in figure 7-10. As mentioned previously, these 
results are not significantly different than the hard bottom results. 

Even after significant trenching has occurred for the 1.25 and 2.0 second excitation pe¬ 
riod cases, the tension results are not significantly different. The evolution of the trenching 
over the first 120 cycles for the 1.25 s case is shown in figure 7-24. The state of the bottom 
after 120 cycles is shown for each of the 1.25, 2.0, and 3.0 s cases in figure 7-25. Corre¬ 
sponding twenty second time series of tension, from the time immediately preceding the 
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Figure 7-24: Changes in the sand bottom over the first 120 cycles of the 1.25 s excitation 
case at At « 0.16. 

(a) T p = 3.0 s, after 120 cycles (b) T p = 2.0 s, after 120 cycles (c) T p = 1.25 s, after 120 cycles 



Figure 7-25: State of the sand bottom after 120 cycles for the (a) 3.0, (b) 2.0, and (c) 1.25 s 
excitation cases at At « 0.16. 


120 cycle mark, are shown in figure 7-26. The presence of the trench increases the mean 
tension level, but it does not change the basic dynamic response. 

The trenching action is a result of the slacking and re-tensioning of the chain following 
a shock discontinuity during the unloading phase of the motion. As the chain re-tensions, 
links on the ground move laterally forward (in the direction of the chain top), carrying sand 
with them. This relationship between the trench and the tension discontinuity explains 
why a trench forms for the 1.25, and 2.0 s cases, but not for the 3.0 s case. In both of the 
higher frequency cases a large pile of sand accumulates at the forward end of the trench. 
In the 3.0 s case, the chain does settle into the sand somewhat, but no pile forms because 
there is no lateral transport of sand by the links. 

These results may have important implications for chain wear in long term deploy- 
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Figure 7-26: Tension time series for the twenty seconds preceding the 120 cycle mark 
on the sand bottom at At ks 0.16, excitation amplitude 0.25 m, and excitation periods 
(a) 3.0 s, (b) 2.0 s, and (c) 1.25 s. 

ments. The lateral motion of the chain along the bottom that is associated with the 
tension shocks may significantly enhance abrasion. If that is the case then wear might be 
reduced by designing moorings so that exceedances of the shock criterion are minimized. 


7.5 Comparison with numerical simulations 

The numerical program described in chapter 3 uses an elastic foundation with linear 
stiffness and damping to model the interaction of the mooring line with the sea floor. 
With the controlled bottom conditions and the importance of the bottom interaction in 
the dynamic response, the laboratory experiments provide an opportunity to investigate 
the limits of the elastic foundation approach in the numerical model. This analysis was not 
possible with the full scale experiment because detailed information about the response 
of the mooring in the touchdown region was not available. 

In searching for a baseline simulation configuration that approximately matched the 
experimental results, it became clear that given correct input for the easily measured 
parameters (mass, weight, static tension), the important parameters in the validation 
were the bottom stiffness and damping, and the chain bending stiffness. Interestingly, 
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none of these three parameters played an important role in the simulations of the full scale 
experiment. The values for these parameters were chosen so that the baseline simulation 
result was in reasonably good agreement with the experimental result. The baseline values 
for these three parameters are k = 10000 N/m 2 , £ = 0.1, and El = 10 -6 Nm 2 . This 
stiffness gives a non-dimensional bottom stiffness, k, of 40.3, 260 times greater than the 
baseline non-dimensional stiffness used for the field experiment. Stiffnesses significantly 
beyond this value made static solutions difficult to obtain. The bottom damping ratio 
chosen was low enough that impact on the tension response is small, but high enough 
that oscillations of the grounded chain are relatively low. The bending stiffness was set 
very low to minimize any possible effect on the response. Hydrodynamic parameters were 
similar to those used for the full scale mooring: Cd n = 0.5, Cd t = 0.01, C an — 0.5, and 
C at = 0-05. 

Figures 7-27 and 7-28 show the simulated response with these parameters for one cycle 
at 1.25 and 3.0 s excitation periods, respectively. In both cases the basic agreement with 
the experimental results on the hard bottom (figures 7-11 and 7-17) is quite good 1 . For 
the high frequency excitation the same snapping and slacking behavior is evident as in the 
experimental result. The magnitude of the tension spike following the snap is higher in 
the simulation, but this may be an artifact of the analog filtering having attenuated the 
impulse in the experimental data. The most significant qualitative difference between the 
simulation and experiment is in the motion of the grounded chain during the downward 
half of the cycle. In the simulation the chain from the rightmost TDP to the current TDP 
at each step is bowed upwards because the model cannot properly resolve the slack in 
the chain along this length. The configuration of the suspended chain (to the left of the 
TDP) does accurately show the reversal in curvature that results from the slack tension 
(or in the simulation, very low tension) in the grounded chain. This discrepancy does not 
appear in the lower frequency excitation case because the tension discontinuity does not 
occur and the grounded chain is not slack. As a result, the simulated profiles in this case 
match the experiment very closely over the entire motion cycle. 

The height of the buckled chain above the bottom during the downward motion can 
be reduced by increasing the damping ratio. Figure 7-29 shows the simulation result for 

1 The experimental results have a small temporal lag relative to the simulation results because of a delay 
in the start of the actuator motion after the video and analog instrumentation is triggered. 
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Figure 7-27: Simulated response with baseline parameters over one cycle at 1.25 s excita¬ 
tion period, 0.25 m excitation amplitude, and At & 0.80. 




Figure 7-28: Simulated response with baseline parameters over one cycle at 3.0 s excitation 
period, 0.25 m excitation amplitude, and At fs 0.80. 
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Figure 7-29: Simulated response with £ = 0.3 over one cycle at 1.25 s excitation period, 
0.25 m excitation amplitude, and At « 0.80. 

C = 0.3 for the 1.25 s excitation case. The motion of the grounded chain is reduced, but 
so are the spatial extent of the motion and the overall tension level. A better alternative 
is to increase the damping ratio to 0.3, but at the same time decrease the bottom stiffness 
to 1000 N/m 2 , so that the damping constant, b = 2C,y/k (m + m a ), remains approximately 
the same as in the baseline configuration. Results for this case, shown in figure 7-30, illus¬ 
trate that lowering the stiffness reduces the height of the buckled chain above the bottom 
while preserving the tension level. This suggests that the bottom damping constant is the 
most important of the bottom parameters in determining the tension and that the motion 
of the chain on the bottom can be largely controlled with stiffness. 

These same parametric variations in the 3.0 s excitation case do not produce significant 
changes in the simulation results. For simulations with (, = 0.3 the maximum tension in 
the cycle changed by a barely detectable 0.08% relative to the baseline simulation. This 
contrasts with the marked decrease in tension and increase in range of motion for the 
1.25 s case. The simulation with £ = 0.3 and k = 1000 N/m 2 at this excitation period 
had similarly small changes. In cases where the shock criterion is never met, the bottom 
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Figure 7-30: Simulated response with ( = 0.3 and k = 1000 N/m 2 over one cycle at 1.25 s 
excitation period, 0.25 m excitation amplitude, and Ar 0.80. 


properties do not appear to play any significant role in the dynamic response. This 
statement comes with the caveat that bottom stiffness always plays a role in the static 
response and thus, to the extent that the dynamic response depends on the steady state 
configuration, the importance of bottom stiffness can never be completely neglected. 

This same situation in which simulation results are much less sensitive to parameter 
variation in the absence of tension shocks is evident in the results with increased bending 
stiffness. With El increased by four orders of magnitude to 0.01 Nm 2 , the results for the 
3.0 s case again only changed very marginally: the maximum tension increased by 1.6%. 
The result for the 1.25 s case is shown in figure 7-31. Both the tension and motion axe 
significantly different than for the baseline simulation. The increased bending stiffness 
allows the wave in the grounded chain to propagate upwards into the suspended chain, 
thus altering the response over the entire length. 

The grounded chain buckles because there is an extended period and region of zero 
tension. With no mechanism to model the collapsing of individual chain links, the chain 
must deform by bending, no matter how low the El value. Providing a means to prop- 
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Figure 7-31: Simulated response with El = 0.01 N/m 2 over one cycle at 1.25 s excitation 
period, 0.25 m excitation amplitude, and At « 0.80. 

agate energy in the presence of zero tension is the very reason for incorporating bending 
stiffness into the equations of motion in the first place. With too high a stiffness, however, 
unrealistic bending effects can propagate into areas with low, but not necessarily zero 
tension. Given the tensions in this model scale system, El = 0.01 N/m 2 is not scaled 
properly to prevent this. This improper scaling is not an issue with the lower frequency 
excitation because flexural waves are never introduced into the system. 

The conclusion of this comparison then is that the elastic foundation is accurate for 
both supersonic and subsonic TDP motions. For the subsonic case this is no surprise given 
the validated accuracy of the simulations of full scale moorings. In the supersonic case, 
there are several qualifiers to this conclusion. Primary among them is the substantial 
sensitivity of the simulation results to parametric variations in bottom properties and 
bending stiffness. This adds additional complexity to the task of defining the simulated 
system. Also, it should be noted that the numerical simulations at the faster excitation 
periods required a higher node density (1601 nodes over the 3.29 m length of the chain) 
than the slowest period cases (401 nodes) to succeed. All of the simulations were run with 
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a 0.001 s base time step. 

Much of the reason for this added difficulty in solutions with tension shocks arises from 
the consequences of the shock condition described in equation 7.11. The non-zero impact 
force gives rise to a dynamic excitation of the elastic bottom. Also, substituting the non¬ 
zero slope at the TDP into equation 7.9 implies a non-zero vertical velocity for the chain 
at the TDP, enhancing the bottom damping forces at that point. To maintain the overall 
accuracy of the simulation, these vibrations must be resolved by increasing the spatial 
and/or temporal resolution of the simulation. The lack of these exciting mechanisms in 
the subsonic case explains the lack of sensitivity of those solutions to variations in bottom 
parameters. 

7.6 Implications for full scale moorings 

With the validation of the elastic foundation approach it is possible to investigate the 
formation of tension discontinuities in full scale moorings, such as that used in the SWEX 
field experiment. As discussed above these tension discontinuities have several implications 
for the design and analysis of these types of moorings. Snap loads and increased wear of 
chain along the bottom may require that design life be shortened or that a heavier material 
be used. 

Tension and TDP and wave speeds for a relatively high resolution numerical simulation 
of the SWEX mooring under the storm conditions of the 3 January 1999 data set is shown 
in figure 7-32. For simplicity, the AxPacks were removed and 401 nodes were used to 
discretize a single 80 m length of chain. The simulation time step was 0.01 s (compared 
to 0.1 s for the simulations in chapter 5 and 6). Snapshots of motion and tension along 
the entire length of the mooring were saved at 0.1 s intervals and used to calculate the 
TDP and wave speeds. Under these extreme conditions the shock condition is exceeded 
during both loading and unloading. The loading shocks correspond to the snap loads 
that are apparent in the time series of top tension (the snap loads are not as clear in the 
experimental tension record because of the analog filtering in the instrumentation). 

To more fully investigate the conditions under which tension discontinuities occur, the 
response of the uniform version of the SWEX mooring was simulated under a range of 
mean tensions and sinusoidal excitations. Current was applied in a linear ramp from top 
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Figure 7-32: (a) Tension at the top of the mooring and (b) TDP speed and transverse 
wave speed at the TDP for a portion of a simulation of the full scale SWEX mooring 
using environmental conditions from the 3 January 1999 storm event. The experimentally 
recorded tension is shown in (a) for reference. 

to bottom, with magnitudes 1.0, 0.6, 0.3, and 0.1 m/s at the top and zero current at the 
bottom. These currents produced non-dimensional mean tensions of 0.245, 0.089, 0.023, 
and 0.003, respectively. Excitation amplitude and period ranged from 0.1 to 2.5 m and 
3.0 to 10.0 seconds. The position and tension at all nodes was recorded every 0.1 seconds 
over the course of each 60 second simulation. This information allows for calculation of 
the position and speed of the touchdown point (in a procedure similar to that used for the 
video data) and a calculation of the wave speed based on the actual instantaneous tension 
at the TDP. 

Figure 7-33 shows the maximum observed difference in the calculated instantaneous 
wave and TDP speeds over the course of the entire simulation during unloading portions 
of the motion. The distinction between loading and unloading motions is made using the 
sign of the TDP speed. Unloading means that chain is being laid down and the TDP 
speed is positive. Exceeding the shock criterion during this portion of the motion implies 
that the chain is being laid down slack. Positive differences in figure 7-33 indicate an 
exceedance of the criterion. The differences are plotted as a function of the ratio between 
the amplitude of the dynamic tension at the top of the mooring and the static tension at 
the static TDP, Tq(0). 
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Figure 7-33: Maximum difference in the wave and TDP speeds during unloading for 
simulations with sinusoidal excitation. At low values of At the results appear clustered 
because there are critical thresholds of input velocity, Au, at which the maximum speed 
difference jumps considerably. 

When the dynamic tension amplitude approaches the static TDP tension, the total 
tension at the TDP approaches zero. With near zero tension the wave speed is very low 
and the shock criterion is easily exceeded. This argument and the results in figure 7-33 
suggest that a reasonable design goal is to keep the value of this ratio below unity to avoid 
tension discontinuities during unloading 2 . As mean tension decreases this goal could be 
relaxed somewhat. Several of the simulation results for the lower values of At are below 
the shock limit but have tension ratios of between two and six. 

The maximum difference between wave and TDP speeds during loading is shown in 
figure 7-34. In this case the difference is plotted as a function of the ratio between the 
input velocity at the top of the mooring (the TDP speed will typically be proportional to 
this) and the wave speed at the TDP calculated from the static tension at the TDP, 

(7.14) 


2 In slack conditions, this goal is nearly impossible to achieve because the static tension at the TDP is 
nearly zero. The effect of the shocks may be slight in these cases, however, because the lack of horizontal 
forces means that the lateral motions at the TDP will be very small. 



TDP and Wave Speed Difference During Unloading 
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TDP and Wave Speed Difference During Loading 



Figure 7-34: Maximum difference in the wave and TDP speeds during loading for simula¬ 
tions with sinusoidal excitation. 

For unloading motion, a tension ratio was used because the shocks are largely dependent 
on low tensions at the TDP. Loading shocks typically form when both wave and TDP 
speeds are non-zero, making a velocity dependence more meaningful. 

As this ratio increases, the TDP is moving faster and faster relative to the wave speed 
and shocks forms. From figure 7-34 an approximate critical value for this velocity ratio is 
0.5. The shock criterion was exceeded for most of the simulations with ratios above this 
value. There are fewer exceedances of the shock criterion during loading than unloading. 
During unloading, 91 of the 140 simulations exceeded the shock criterion (105 simulations 
have a tension ratio greater than unity). During loading, only 41 simulations exceeded 
the shock criterion (47 simulations have a velocity ratio greater than 0.5). That loading 
shocks (snap loads) are more rare than unloading shocks (slacks and lateral motion along 
the bottom) is consistent with the experimental results where the results for 3.0 second 
excitation period had no shocks, 2.0 second period had unloading shocks, and only the 
fastest excitation cases had both loading and unloading shocks. No simulations in the 
above cases had only loading shocks. 

The critical values of the tension and velocity ratios described above can be used along 
with the simple model for dynamic tension described in chapter 6 to develop a procedure 
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for estimating the likelihood of tension discontinuities in full scale moorings. Given an 
input wave spectrum, spectra of heave velocity and acceleration can be computed and used 
in equation 6.31 to calculate a spectrum of tension at the top of the mooring. Ass um ing 
that the tension is a Gaussian random process, the expected number of times that the 
dynamic tension will exceed the TDP static tension, T 0 (0), per second is [73] 


ML e -TZ(0)/2Mf 

2tv y Mq e 


(7.15) 


where Mg' and Mj are the moments of the tension spectrum, 


poo 

Mq — / Sx(uj)d(jj — <t|,, 

Jo 

(7.16) 

M2 = / LJ 2 ST{w)dtV. 

Jo 

(7.17) 


The probability of at least one exceedance in a period t is 


P(|T| > T 0 (0) in t) = 1 - e~ iNT . 


(7.18) 


Similarly, for loading shocks the number of exceedances of the input velocity of the level 
|lwave(0) per second is 


7\T — _^. e - y wave(°)/ 8 ^? 

ly v — 0 e i 

cr v 

and the probability of at least one exceedance in t is 


(7.19) 


P(Auj > 0.5V wave (0) in t) = 1 — e tNv . 


(7.20) 


In the above, the variances of velocity and acceleration have been substituted for the zero- 
and second-order moments of velocity. 

To test these guidelines, this same mooring configuration was simulated with the steady 
state and dynamic excitation conditions observed from the SWEX experiment. Each of 
the 119 simulations was run for 200 seconds (t — 200 s). The moments for tension 
were calculated using equation 6.31 with model coefficients calculated from a fit to the 
simulation results (M = 158.0 kg, Cd = 0.288). From the simulation results, 10 data 
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shock type 

correct + 

correct — 

false + 

false — 

unloading 

112 

0 

7 

0 

loading 

4 

109 

0 

6 


Table 7.2: Number of correct and incorrect predictions given a probability level of 0.9 in 
equations 7.18 and 7.20 as an indicator of the presence of shocks. 

sets had loading shocks and 112 simulations had unloading shocks. Using a probability 
level greater than 0.9 as an indicator that a shock could be expected, the accuracy of 
equations 7.18 and 7.20 can be categorized in one of four ways: correct positive, correct 
negative, false positive, and false negative, where a positive is an exceedance of the shock 
condition. For example, a correct positive is a data set for which the calculated probability 
of an exceedance is greater than 0.9 and an exceedance was observed in that data set. 
Table 7.2 lists the number of data sets that fall into each category for both types of shock. 

There are no negative predictions for the unloading case. The calculated probability 
of the tension ratio exceeding unity is nearly 1.0 for all data sets. The seven data sets 
in the simulation results that did not have an exceedance during unloading were amongst 
the twelve lowest of all data sets ranked in terms of dynamic tension (cr T ). Unloading 
shocks occurred under nearly all of the observed conditions. In this situation then, the 
probabilistic prediction of unloading shocks in all 119 cases is not unreasonable. 

As observed in the study with sinusoidal inputs, loading shocks occur less frequently 
than unloading shocks. Equation 7.20 with a probability level of 0.9 appears to offer 
a reasonable predictive capability for loading shocks. However, for conservative design, 
the false negatives are a concern. The number of false negatives can be reduced, with a 
subsequent increase in the number of false positives, by decreasing the probability level. 
A value of 0.5 produces only two false negatives, but also yields three false positives. 

Overall then, the design guidelines outlined above provide a reasonable prediction of 
the likelihood of tension discontinuities at the TDP. For the most accurate prediction, and 
for quantitative information about the magnitude of tension spikes and extent of lateral 
motion along the bottom associated with shocks, full numerical simulation (with accurate 
representation of bottom conditions) is still necessary. 
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Chapter 8 


Conclusions 


The most tangible contributions of this thesis are tools that can be used in the analysis and 
design of mooring systems. The generalized-a time integration scheme and algorithms for 
mesh refinement, adaptive time-stepping, and adaptive relaxation contribute to the nu¬ 
merical program and make it robust and relatively easy to use. The simple model for 
dynamic tension in chapter 6 can provide a mooring designer with a convenient and accu¬ 
rate predictor of tension given very simple inputs. On a more fundamental level, however, 
the tools are not themselves the end goal of this work. That goal is to develop a deeper 
understanding of the mechanics of these systems so that design methodologies can be im¬ 
proved, and more capable, longer lasting systems can be developed and deployed. Toward 
that end, the real importance of the tools is the insight that they yield in the analyses 
such as those of dynamic tension in chapter 6 and bottom interaction in chapter 7. Tools 
facilitate design and analysis, but ultimately, innovation must come from understanding. 

8.1 Summary 

8.1.1 Numerical model 

The generalized-a time integration scheme for cable dynamics developed in chapter 2 
offers significant advantages over the traditional box method [1] or other box method 
variants [60]. By retaining the box method’s finite difference spatial integration, the 
method remains second-order accurate in the spatial dimension and is relatively easy to 
implement. For the temporal discretization the generalized-a algorithm provides: 
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• Controllable numerical dissipation without loss of second-order accuracy. Trape¬ 
zoidal rule is only first-order accurate in the presence of numerical dissipation. The 
original box method temporal discretization has no numerical dissipation and there¬ 
fore is subject to Crank-Nicolson noise and other numerical instabilities. 

• Second-order accuracy. Box method variants using backward differences are only 
first-order accurate in time. 

• The ability to implement other algorithms, including backward differences, trape¬ 
zoidal rule, HHT-a, and WBZ-a, through appropriate choices of parameters. 

• Improved numerical stability through the averaging of coefficient matrices. 

These advantages were made clear in chapter 5 where the scheme allowed for robust 
solution of the instability in the two-dimensional hanging chain motion leading into three- 
dimensional whirling. This solution could not be obtained with the previous program 
that used the pure box method. The generalized-a scheme also facilitated the fast, accu¬ 
rate, and robust simulation of the entire range of conditions observed during the SWEX 
experiment in chapters 5 and 6. 

8.1.2 Models for understanding dynamic tension 

Chapter 6 contains a number of significant contributions regarding dynamic tension in 
geometrically compliant systems. The proposed model for dynamic tension was derived 
by first considering a SDOF spring-mass-dashpot system. By fitting motion and tension 
spectra for each data set individually to this form, model coefficients were derived for 
each data set. To understand the scatter in these coefficients, the individual tension 
mechanisms (inertia, drag, stiffness) were analyzed separately and in pairs. This analysis 
confirmed that the model coefficients should change roughly linearly with mooring shape 
(as measured by the non-dimensional mean tension) and that the different mechanisms 
are coupled together. 

To account for scattering due to this coupling, relationships were sought among the 
statistics of the experimental data that presented low scatter. From these low scatter rela¬ 
tionships the tension model was constructed based on the standard deviations of tension, 
heave acceleration, and heave quadratic velocity. Analyzing the model in relationship to 
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the spring-mass-dashpot model showed that coupling between mass and drag was a likely 
cause of much of the scatter in the coefficients fitted to individual data sets and that the 
spring terms could be neglected over the range of conditions present in the experimental 
data set. With just two parameters the simple model is compact and able to represent the 
experimental data with only slightly less accuracy than the model using the coefficients 
fitted to individual data sets. 

Numerical simulations were then used to analyze the parametric dependence of the 
model mass and drag coefficients on actual system parameters, including normal and 
tangential drag coefficients, added mass coefficients, and bottom stiffness and damping. 
From these parametric dependencies, formulae were derived for the a priori prediction of 
the model coefficients. These formulae, together with the simple model, allow a designer 
to analytically calculate the dynamic tension response over a wide range of conditions 
without the need for experimental data or numerical simulation. This approach was 
validated using data from a second oceanographic mooring (the CMO experiment) and 
simulation of steel catenary and lazy wave riser configurations. 

Several circumstances under which the simple model dynamic tension is not accurate 
were described. Among these are: 

• The coupling between mass and drag is self-limiting. At high sea state, the coupling 
term in the model leads to an over prediction of the dynamic tension. 

• At very low frequencies the only dynamic tension is due to stiffness effects which are 
neglected in the model. These effects are small. 

• There is no model drag when the mooring is slack/vertical. If input velocities are 
high in such a case, the dynamic tension calculated from the model will be too low. 

• The model is derived based on an assumption that the mooring materials are inex- 
tensible. If the mooring has significant elastic compliance then the model calculated 
tensions will be too high. 

• The model assumes that adequate scope is always available for geometric compliance. 
When the mooring is pulled taut, model calculated tensions will be too low. 

An additional limitation of the simple model is that as Ar, the non-dimensional mean 
tension, increases, dynamic tension due to horizontal motions becomes significant. In the 
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SWEX mooring, the non-dimensional mean tensions were low enough that surge and sway 
motion did not contribute significantly to the total dynamic tension. For oceanographic 
moorings in very shallow water, however, horizontal motions become important, as illus¬ 
trated by the data from the NDBC Duck Pier mooring in 17 m depth. This situation 
is also more common in riser applications where mooring pre-tensioning can lead to very 
high values of At. 

In the statistical domain, the contributions to the dynamic tension due to vertical 
and horizontal motions are linearly separable. Chapter 6 demonstrated that a model for 
dynamic tension due to horizontal motion can be developed using a procedure similar to 
that used for vertical motion. Summing the results from the two models yields a complete 
prediction for the dynamic tension in the presence of both vertical and horizontal input 
motions. 

8.1.3 Bottom interaction 

The experiments presented in chapter 7 represent the first direct observation of the shock 
condition at the touchdown point of a catenary mooring. They also illustrate the impli¬ 
cations of the shocks. The mathematical implications of the shock are a non-zero impact 
force at the TDP and a loss of tangency as the mooring line leaves the bottom. Practi¬ 
cally, the shock condition has different implications depending on whether the motion is 
loading, line being picked up off the bottom, or unloading, chain being laid down on the 
bottom. 

In the unloading portion of the cycle, the tension discontinuity leads to zero tension 
at the TDP and the grounded line being laid down with slack. As the motion reverses the 
mooring line does not roll smoothly off the bottom but rather the TDP moves laterally 
along the bottom. This lateral motion may be a significant source of abrasion of mooring 
line in the touchdown region. This type of tension discontinuity is more common than 
shocks during loading. 

During loading, the tension discontinuity of the shock is manifest as a snap load, leading 
to large impulsive tension spikes. This situation arises only after the shock criterion has 
been exceeded during the unloading motion. The snap occurs because the geometry of 
the mooring cannot change rapidly enough to accommodate the retensioning of the slack 
grounded line. 
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The qualitative features of the response of the mooring in the presence of tension shocks 
at the TDP do not change as the bottom type changes. On two softer artificial bottoms 
with higher friction, there were virtually no observable quantitative or qualitative changes 
in the motion or tension response. On a realistic sand bottom, the lateral motion of the 
TDP associated with unloading shocks caused a large trench to form in the touchdown 
region. The digging action required to produce this trench strongly supports the idea that 
mooring wear and abrasion will increase in the presence of unloading shocks. 

The likelihood of both types of shocks can be predicted using relatively simple criteria. 
However, for the most accurate prediction of their occurrence and their implications, 
full numerical simulations should be used. The presence of the shocks make numerical 
simulations more difficult because the results are more sensitive to parametric variations 
in the bottom parameters and the mooring bending stiffness. The results from chapter 7 
did verify that with the right parameters, the elastic foundation approach to modeling 
bottom interaction can capture the tension and motion response in the touchdown region 
quite accurately. 

8.2 Recommendations for future work 

In its current state of development the numerical program is relatively robust and capable. 
While not a focus of this thesis, additional work is still needed to improve the stability 
of static solutions, particularly for problems with mooring line on the sea floor. The 
current static solution procedures are adequate, but a fast, robust, fully automated scheme, 
coupled with the generalized-a based time integration scheme developed in this thesis 
would be a very powerful tool for mooring line simulation problems. 

The instrumentation suite that was deployed for the field experiments provided good 
quality data for the validation of the numerical program and for the analysis of dynamic 
tension. It did not provide data that could be used for a full scale comparison with 
the results from the laboratory experiment. For that analysis, more complete information 
about the steady state configuration and along mooring motion are needed. Redeployment 
of the GPS receiver that failed in the SWEX experiment and an acoustically localized 
anchor position would yield high quality data about the very slow current and tidal induced 
motions of the mooring. The AxPack instruments would provide much more information 
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about motion over the length of the mooring with the addition of a compass and tilt sensor. 
In conjunction with the accelerometer data, this addition would provide high frequency 
earth referenced motion along the mooring. 

For the laboratory experiments additional realistic bottom types need to be studied. 
More work is also needed to refer the laboratory observations to full scale bottom condi¬ 
tions. Bottoms such as mud will likely make the video data very difficult to process and 
additional analog instrumentation, such as very small inline accelerometers and load cells, 
may be necessary to develop an accurate picture of the system dynamics on these bottoms. 
An impact resistant, unobtrusive, inline load cell would also be useful in collecting tension 
data directly in the grounded portion of the chain. 

An interesting laboratory experiment would also be one in which chain abrasion in 
the presence of tension shocks could be measured and compared to the wear experienced 
during more typical, smooth motions. Such an effort would require considerable thought 
about the practical aspects of runs that may last for many days or weeks. Together with 
design tools that can predict the occurrence of tension shocks, a catalog of this kind of 
data on various bottoms would be very valuable, as chain wear is one of the limiting factors 
in current oceanographic mooring practice. 

As discussed in the summary for chapter 6, additional work is needed to develop 
design formulae for dynamic tension in compliant systems with horizontal input motions, 
analogous to equation 6.16 for vertical motions. A reduction to such a simple form may 
not be possible, but through the study of a range of systems in which horizontal motions 
are important some general design rules could certainly be formulated. 
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Appendix A 


Derivation of 2D Equations of 
Motion 


This appendix contains a derivation of the two-dimensional governing equations for a cable 
in water. A derivation of the three-dimensional equations can be found in Tjavaras [93]. 
The derivation assumes that the cable material is circular and homogeneous in cross- 
section (but not necessarily along the length), has a nonlinear tension-strain relationship 
and that Euler-Bernoulli beam theory can be applied. Fluid forces on the cable are 
modeled using a Morison formulation [29]. 


A.l Kinematics and coordinate system 


The governing equations are derived in the coordinate system defined by the local tan¬ 
gential ( t) and normal (n) directions, as shown in figure A.l. The transformation between 
local and global (i, j) coordinates is 
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Figure A-l: Vector definitions for the local coordinate system. 
The time derivative of a vector, A, that is defined in the local frame is 


dA 

dt 


dA 

dt 


“b D X 


A, 


(A.3) 


where u) is the time rate of change of the orientation of the local frame. Likewise, the 
derivative of A with respect to the Lagrangian coordinate, s, along the cable length is 


dA 

ds 


dA 

ds 


+ fi X 


A, 


(A.4) 


where f2 is the spatial rate of change of the orientation of the local frame. For the two- 
dimensional case defined in figure A.l, u) and 0, are 

a -£s- ^=T k - (A - 5) 


A.2 Cable stretch and buoyancy 

If ds is the unstretched length of an infinitesimal element of cable and dsi is the stretched 
length then 


dsi = (1 + e)ds, 


(A.6) 
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where e is the cable strain. Prom conservation of mass, the mass and weight per unit 
length of the stretched element are 


mids\ = mds —> m\ = 
W\ds\ = wds —»• tni = 


m 

TT? 

w 

1 + € 


(A.7) 

(A.8) 


Based on Poisson’s ratio, v , the reduction from the nominal diameter, d, of the stretched 
cross-section is 


Sd = (—ve)d, 


(A.9) 


the change in cross-sectional area is 

SA = ^ [(d + Sd) 2 - d 2 ] ss ~dSd, (A. 10) 

and thus 

SA 

= “ 2z/e - (A-il) 

If v = | we have conservation of volume 1 

A{ 1 — e)(l + e)ds ~ Ads. (A.12) 

Finally, we can use a binomial expansion to write the stretched area and diameter in a 
more convenient form: 

M = A(1 - e) « (A.13) 

iff 


1 While v = | is not strictly true for all cables (particularly metal chain and wire), the conservation 
of volume that it introduces greatly simplifies the treatment of the buoyancy forces on the cable [36]. 
Burgess [10] calculates the error associated with using v — | when the true value of v = 0 (the 
maximum possible error in u) as pgz/E. This term becomes significant only with large depths and very 
soft materials. 
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Reality is pressure on Reality + pressure applied Compensating tension vectors 

element exterior only at the two ends 


Figure A-2: Schematic diagram of pressure and effective tension terms. 


and 


di 


d 

vm 


(A.14) 


With the above definitions, we can easily treat hydrostatic forces on the cable by 
considering the effective tension. In reality, hydrostatic forces act only on the exterior of 
the element, not at the two ends. Following the procedure first suggested by Breslin [6], 
however, we can introduce a pressure force on the element end faces if we also introduce 
a compensating term into the tension force. This is shown schematically in figure A-2. 
Mathematically, 


^effective T pA \, 


(A.15) 


where p is the hydrostatic pressure at the depth of the element. With the fictitious end 
pressures Archimedes’ principle applies and the buoyancy force per unit length of the 
stretched element is simply 

Bi = Aip w g = ——Pw9, (A.16) 

i ~r € 


where p w is the density of water and g is the local acceleration due to gravity. The total 
of the weight and buoyancy forces on the stretched element are 



(l + e Pw9 



(A.17) 
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If we define the wet weight of the material as wq = w — Ap w g, then 


(1 + e)F w = — wqi = — wq cos fit + wq sin cpn. (A.18) 

Introducing the effective tension and the wet weight frees us from any further con¬ 
sideration of the hydrostatic pressure; pressure effects are now simply rolled into any 
computed strain result [36]. Because tensions are always computed as a function of strain, 
all calculated tensions will be the effective tension. For simplicity in the remainder of the 
derivation of the governing equations we will use T to denote the effective tension. 

A.3 Hydrodynamic forces 

The hydrodynamic forces on the cable are the drag, added mass, and dynamic Archimedes 
forces. The drag forces arise from the relative velocity of the cable in a current field 
defined in global i, j coordinates by U and V respectively. In local coordinates the relative 
velocities are 


u r — u — U cos (f) — V sin <f>, and (A.19) 

v r = v + U sin (j> — V cos 4>. (A.20) 


The drag force in local coordinates is 


F d = 


2 P w yT+e ^dt ( u u c) | U U c \ t 

~ 2P w i/E+e ^d„ ( v ~ v c) \ v ~ v c\ n 


(A.21) 


For a solid circular cross-section cable, the added mass force has a component in the 
normal direction only. It is computed as a function of the relative acceleration between 
the fluid and the cable. The time derivative of the current velocity in local coordinates 
(assuming steady current) is 


v c 


— [U cos <p+V sin (j)\ 



(A.22) 
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The added mass force is 


m„ 


1 + e 


d(f> dv 

(U cos cf) + V sin <f>) — - — 


n. 


(A.23) 


where m a is the added mass per length of the cable cross-section. Because the current 
velocity in local coordinates is changing in time, there is a pressure gradient that gives 
rise to the dynamic Archimedes force [69]. Like the added mass force, the only component 
of this force on a solid circular cable is in the normal direction. It is defined as 


Far = -pw ^ (U COS 0 + V sin <f>] ‘—n. (A.24) 

For cables with non-solid, or non-circular cross-sections (such as chains) there can be 
both added mass and dynamic Archimedes forces in the tangential direction. The time 
derivative of the current velocity in the tangential direction is 

ii c = [-U sm4>+V cos (A.25) 

With the appropriate tangential components equations A.23 and A.24 become 


F a m. — 


1+e 

m an 

L 14"C 


\—U sin <f> + V cos <£) — fijr 

— (U cos (f> + V sin 4>) §* — ff 


f? mg -wo 
ar ~ g(i + e) 


(-U sin <t> + V cos (j)} | 

— (U cos (f> + V sin <j)} ^ n . 


(A.26) 

(A.27) 


The term defines the mass of the fluid displaced by the irregular cross-section. 

This formulation also requires two terms to describe the cross-section added mass, m at 
for tangential motion, and m a „ for normal motion. 


A.4 Balance of forces 


A summation of the forces on the cable element shown in figure A.l yields 


A 

dt 




= T + df-T+Fds 1 . 


(A.28) 
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If we expand the derivatives according to equations A.3 and A.4 and eliminate stretched 
variables we find 


dT 

m - 7 ——|- w x V - ——|- x T + ,F(1 + e). 


(A.29) 


Substituting F = F am + F a d + F w + Fa and collecting terms in the normal and tangential 
directions yields 
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(A.31) 


A.5 Balance of moments 

For the two-dimensional element, the only moment balance involves moments about the 
out-of-plane axis. The tension force does not contribute a moment in this case because in 
the infinitesimal limit the tangential and normal directions at the two opposite ends have 
the same direction and opposite magnitudes. The remaining moment contributions are 
the rotational inertia, the couple due to shear, and the bending of the element: 


d ( dd>\ 

— \dsipJi— 1 = S n ds 1 + dM u 


(A.32) 


where I\ is the second-area moment of inertia of the stretched cable cross-section and p c 
is the mass density of the cable. Using Euler-Bernoulli beam theory the bending moment, 
M, is the product of the flexural stiffness of the cable, El, and cable curvature: 


Mi = Eh- 


(A.33) 


The area moment of inertia of a circular cross section is ~ d 4 and thus 


M = Mi (1 + e ) 2 7 = J 1 (l + e) 2 . 


(A.34) 
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Substituting equations A.33 and A.34 into equation A.32, eliminating ds i, and dividing 
by ds yields 


£ ( Pd d<A 

dt \ 1 + e dt ) 


_5 n (1 + e) + £ 
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\(1 + e) 2 ds ) ' 


(A.35) 


Expanding derivatives and re-inserting the definitions for the spatial and temporal deriva¬ 
tives of (j) given in equation A.5 so that the system remains a first-order PDE yields, 


Pci 


(1 + e) 


Ouj de 


dt 


iv 


dt 


= El 


dCl 3 

ds 


2fl 3 de 
1 + e ds 


+ (1 + e) 3 S n . 


(A.36) 


Howell [46] used dimensional analysis to show that for both metal and synthetic cables 
the rotational inertia term and the second bending term containing the spatial derivative 
of strain are of significantly lower order than the remaining terms. If we drop these terms, 
the result is 


El 


dfl 3 

ds 


+ (1 + e) 3 S n . 


(A.37) 


Both terms could be retained without a significant loss of simplicity. Using equation A.37 
over equation A.36 does offer the advantage that the additional dependent variable u does 
not need to be stored. 


A.6 Compatibility 


Compatibility can be established by requiring continuity of the position of the cable in 
both space and time. If R(s, t ) is a vector to a point on the cable then continuity requires 


£ 

dt 



(A.38) 


By definition 



(A.39) 


and from analytic geometry we know that the derivative of a position vector to a space 
curve with respect to arc length is a unit vector tangent to the curve, in the direction of 
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increasing arc length [43] 


dR ~ dR - 

d^ = t ^di - (1+e) ‘- < A ' 40 > 

Thus the continuity condition (equation A.38) reduces to 

Expanding the derivatives and collecting components in the tangential and normal direc¬ 
tions gives 


du de dip 
ds dt V d s' 


dv 

ds 


(1 + 0 


d(j> 

dt 



(A.42) 
(A.43) 


A.7 Matrix form of the governing equations 

Equations A.5, A.30, A.31, A.36 or A.37, A.42, and A.43, define a system of either six 
equations and six unknowns (without rotational inertia) or seven equations and seven 
unknowns (if rotational inertia is retained in equation A.36). The six degree-of-freedom 
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form of the equations can be written as 



' n ds 


du dp 

— m— + mv— - wo cos p 

dt at 


^PwdnCdtUr |uj-| 's/l T c 0, 


(A.44) 


dS n 

ds 


+ T(e)^ - (m + m a ) 
as 


dv 

dt 


mu + 



(U cos cf) + V sin p) 


d(j) 

dt 


+ w 0 sin p - -p w dC dn v r |v r | VTTe = 0, 


(A.45) 


du d(f> de 

Ys~ V d~s~di^ ’ 

(A.46) 

dv dcj) .d<p 

9j + “8J-< 1 + £ > Yt=°’ 

(A.47) 


(A.48) 

El ^ + S n (1 + ef = 0. 

(A.49) 


If we define Y = [e, S n ,u,v,p, fIs] T then equations A.44 through A.49 can be written in 
matrix form as 


M ^ +K ^ + F(Y ' s ' t) = 0 ' 


(A.50) 


The continuous forms of the mass matrix, “stiffness” matrix, and forcing vector are 


M = 


0 0 —m 

0 


0 0 

-1 0 

0 0 

0 0 

0 0 


0 

0 

0 

0 


0 

~{m + m a )j -1 
0 
0 
0 
0 


mv 

| mu + ( p w ^ + m aj (U cos <f) + V sin p) 
0 

-(1 + e) 

0 

0 


(A.51) 
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0 


K = 


T’(e) 0 0 0 S n 

0 10 0 T(c) 

0 0 1 0 -v 

0 0 0 1 u 

0 0 0 0 1 

0 0 0 0 0 


0 

0 

0 

0 

El 


(A. 52 ) 


F = 


—wq cos (f) — | p w ndCd t u T |u r | %/l + e 
w 0 sin <f>+ \p w dC dn v T |u r | %/l + e 


0 

0 

—Q3 


<Sn (1 + e) 3 


(A. 53 ) 


Note that the distribution of terms as either stiffness or force is somewhat arbitrary as 
spatial derivatives of <j> could also appear as O3. Experience has shown that given the 
dependence of the nonlinear solver (described in appendix C) on the Jacobian of the 
equations of motion, the solution typically proceeds more quickly when any dependence 
on (j> is explicitly incorporated into the equations of motion. This is not surprising given 
the formulation outlined above in which cf> is the primary variable used in describing the 
system geometry. 
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A.8 Static governing equations 


The static governing equations can be derived from the dynamic equations (A.44 to A.49) 
simply by dropping time derivative and velocity terms. They are 


T'(e)~ - S n ~ - wo cos <f)+ ~p w dnCd t (U cos cj) + V sin (f>) \U cos (f) + V sin0| VTTe = 0, 
as os 2 

(A.54) 


+ T(e)^ + wo sin (f)+ \p w dCd n (-U sin <f) + V cos <f>) \—U sin </>+V cos</>| Vl + e = 0, 
os os 2 



(A.55) 

— O 3 = 0 , 

(A.56) 

^ + S n (l + e ) 3 = 0 . 

(A.57) 


226 


Appendix B 


Accuracy and von Neumann 
Stability Analysis of the Box 
Method 


B.l Stability 

In this appendix we use the classical von Neumann method [48,82] to analyze the stability 
of the box method as a pure finite difference method. This contrasts with the amplification 
matrix method of stability analysis which operates on the semi-discrete equation of motion. 

Like the amplification matrix method, we consider a single degree-of-freedom homo¬ 
geneous problem 


dy 

dt 


dy 

+ “<l= 0 ' 


(B.l) 


The fully discrete form of this equation after applying the box method is 


Vi - y 


n~ 1 


At 


+ 


_ ,, n ~ 1 

Vj-i yj -1 

At 


1 

+ 2 W 


qiV" — 77 /^ ^ _ 7/^ ^ 

y j yj -1 , y± y j -1 

As As 


= 0 . 


(B.2) 


The proposed solution is 


y n = c n e ijkAs ^ 


(B.3) 
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where ( is an amplification factor, i = \/—1, and A; is a spatial wave number. The condition 
for stability of the method is |£| < 1. 

Substituting equation B.3 into equation B.2 and dividing through by ( n ~ 1 e ^ kAs yields 

= (l+e-' ta -)-u4f (l-e-* a ‘) 

^ (1 + e ~ ikAs ) + u;ff (1 - e - ikAs ) ' 

For all values of k, As, At, and u, |(j = 1, and as in chapter 2 we find that the box 
method is unconditionally stable with no numerical dissipation. 


B.2 Accuracy 


The truncation error of the box method is found using a procedure similar to that for the 
semi-discrete equations in chapter 2. Given an exact solution to equation B.l, y, Taylor 
series expansions for the solution near s = jAs and t = nAt are written as 


vY 


~n 

Vi -1 


~n—1 

Vj-1 
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= yj-As 


dy\ n 4. (&vY _ (&v\ n 
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+ ... 
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dtjj \ ds / j 
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At 3 /d 3 y\ n As 3 / d 3 y\ n At 2 As f d 3 y \ n AtAs 2 

w), ~ IT yd^Jj 2 2 


Because y is an exact solution to equation B.l we can write 


a S) n = -J a A 

ds) j \ds, j 

&y\ n -- (*LY 

dfijj U \dsdt)j ’ 


( d 2 y V 
\ dsdt) 

( d 3 y \ n 
\dsdt 2 J j 


= — U3 
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( d s y V 

W \ds 2 dt) j ' 


(B.5) 

(B.6) 

(B.7) 


( d s y Y 

\ds 2 dt) • 


+ ... 


(B.8) 

(B.9) 

(B.10) 

(B.ll) 
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Using these relationships and substituting the Taylor expansions for the exact solution into 
the approximate difference equation B.2 yields an expression for the truncation error [83], 


1 

2 

At 2 


f," _ _ f.n-l 

y j y_j + 9j -1 2/j-l 


At 

1 /£yV 


At 
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f/ n — f/ n ri 71-1 
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6 ydfi)- + 4 \dsdt 2 ) ■ 
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( ° 3 y v 


+ As 2 


() 


+ 


a> 


4 \ds 2 dt J • 6 \ ds 3 ) ■ 


sW 


+ H.O.T. (B.12) 


The truncation error that results from using the difference equation in place of the exact 
PDE has terms of lowest order in At 2 and As 2 . Thus, the method is second-order accurate 
in both space and time. 
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Appendix C 


Solution of the Nonlinear Problem 

For all combinations of boundary conditions, 2D or 3D and static or dynamic problems, 
the mathematical problem is posed as a system of coupled, nonlinear partial differential 
equations. These systems are solved numerically by discretizing the continuous (exact) 
forms of these governing equations onto a grid of nodes and calculating an approximate 
solution. As the grid becomes finer and finer the approximate solution will approach the 
exact solution. The cost of these finer discretizations which buy better solutions is an 
increase in computation time. 

Both the static and the dynamic cable problems can be generalized as a system of N 
first-order nonlinear partial differential equations (at each time step the dynamic problem 
represents a quasi-static equilibrium problem), 

K^ + F(s,Y) = 0, (C.l) 

where Y is the vector of the N dependent variables. For example, in the 2D static prob¬ 
lem (the simplest of all possible cases), equation C.l represents four equations in four 
unknowns: strain (from which we can always derive tension via a constitutive relation¬ 
ship), shear force, inclination angle, and curvature. This equation is discretized at the n 
nodal points using centered finite differences written on the half-grid points. At node j 
the discretized result is 

Yj — Yj-i + —— ~~~~ (F j + Fj_x) = 0. (C.2) 


231 



When combined with a total of N boundary conditions at the two ends, equation C.2 
written at the n — 1 half-grid points represents a coupled system of N x n nonlinear 
equations in N x n unknowns. The system can be solved using a relaxation procedure 
similar to Newton-Raphson [82]. 


C.l Newton-Raphson updates 

Equation C.2 can only strictly be satisfied by an exact solution for Yj. Given an inex¬ 
act first guess at this solution, Y°, we need to develop an iterative scheme to calculate 
successively better approximations, Y*, through a series of update vectors, AY*-, such that 

Y i+i = Y *. + AY*-, (C.3) 


where Yj +1 brings us closer to satisfying the equality in equation C.2. In quantitative 
terms we want to iteratively minimize the error function 


eJfYi.Y}-,) =Y]-Y5_j + 


Sj Sj-l 


(1 + Fj-i) • 


(C.4) 


Neglecting for clarity the dependence on the previous nodal point (j — 1), we can very 
loosely write 

ey(Y; + AY-)- e ‘(Y-)] _ a,. 

AY*. ~ 8Yj' 

The derivatives on the right hand side of equation C.5 can be calculated analytically from 
the known form of the discretized governing equations (equation C.4). If we were to re¬ 
insert the dependence on Yj-i, we would note that these derivatives actually constitute 
an N x 2 N Jacobian matrix at each j (the matrix is composed of the derivatives of the N 
equations with respect to the 2 N variables represented by Yj and Yj_i). We can assemble 
the Jacobian matrices from each node into a single global Jacobian matrix (much like 
element stiffness matrices are assembled into global stiffness matrices in the finite element 
method), add boundary condition information and formulate a linear system that will 
find AY*- to drive the updated error, e* +1 , to zero. If J* is this global Jacobian matrix 

3 J 
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evaluated at Y 1 then we see from equation C.5 that 

J'AY* = -eh (C.6) 

Because only two nodes (j and j — 1) are coupled by each individual Jacobian matrix 
the assembled global Jacobian matrix in equation C.6 will be very sparse, with the only 
non-zero entries clustered near the main diagonal. We can take advantage of this sparsity 
in solving equation C.6 by using a sparse Gaussian elimination algorithm, NSPIV, due to 
Sherman [85]. Sparse algorithms such as NSPIV exploit sparsity to reduce both memory 
requirements and computation time (normal Gaussian elimination is an 0(n 3 ) operation, 
sparse algorithms can be as efficient as O(n)). A distinct advantage of NSPIV over some 
other algorithms for sparse linear systems is that it can handle matrices with arbitrary 
sparsity patterns. This capability is important in dealing with systems that are not simply 
connected (i.e., multipoint moorings and moorings with segments that branch out from 
other segments). 

C.2 Convergence criterion 

The iterative updates of Y continue until the updates, AY, become sufficiently small as 
to not warrant continuation of the process. The total error at iteration i, cr *, is defined as 



where AY j k are the N components of AY at node j and iteration i and Xk are scaling 
constants appropriate to each of the physical variables represented within Y. The stopping 
criterion is simply 


a 1 < specified tolerance. 


(C.8) 
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C.3 Relaxation 


The actual update to Y is scaled by a relaxation factor /i r 

Y i+l = Y* + Hr AY*. 


(C.9) 


The purpose of this relaxation factor is to slow (under-relax) the update in cases where 
strong nonlinearities may mean that the update is not quite as robust as we would like. 
For highly nonlinear problems, where small changes in parameters can mean large changes 
in system configuration, the approximation of equation C.5 becomes less valid and our 
update AY*, if fully applied (p r = 1), may actually increase the total system error. A 
small relaxation factor increases the accuracy of the linearized Taylor series expansion 
that equation C.5 represents. By slowing the process down (p r < 1) the movement of the 
system from iteration to iteration towards equilibrium will be smoother because the steps 
between iterations will be smaller. 

In many cases, it is desirable to have the relaxation factor vary as the solution pro¬ 
gresses. This is particularly true in the static solution of some problems which may need 
very fine movement of the relaxation process as the solution approaches equilibrium. As 
an example, in cases with cable resting on the sea floor, the resolution of the location of the 
touchdown point can be very difficult because of the unilateral nonlinearity represented 
by the bottom. With too large a relaxation factor the update might pull a substantial 
amount of cable off the bottom, only to be followed by an update that drops too much 
onto the bottom. 

To accommodate this behavior the actual relaxation factor used from iteration to 
iteration is varied according to the progress of the solution. If at any point during the 
solution the error increases from one step to the next, a z > cr* _1 , the relaxation factor 
applied to the update at that step is reduced from the factor used at the previous step, 


^ -l- 

(7 > a 


Mr 


/4 1 

Ri : 


(C.10) 


where R\ is a constant larger than unity. If the error is decreasing as it should then we 
take the opposite approach and try to speed the solution by increasing the relaxation 
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factor 


u 2 ' < cr* 1 -> 4 = fia/i 4 " 1 , (C.ll) 

where 1 < i ?2 < -Ri • The relaxation factor is not allowed to increase beyond the baseline 
value, /x r , and as a protection against pathological cases in which a very small relaxation 
can effectively bring the solution to a halt, it is not allowed to decrease beyond fj, r / 1000 . 
In our implementation Ri =1.1 and R <2 = 1.02. 

The adaptive procedure has the effect of driving the relaxation factor into an equi¬ 
librium at which the solution can make the best progress. For regions of the solution in 
which the baseline relaxation is too large and the solution starts to diverge, equation C.10 
kicks in and the relaxation factor is reduced until the solution begins to converge again. 
Equation C.ll mitigates these reductions and prevents the relaxation factor from getting 
too small in response to an occasional wayward oscillation in an otherwise downwards 
solution path. 

This procedure still requires a reasonable value for the baseline relaxation factor, but it 
avoids having to set that factor very low when in fact it may need to be very low only for a 
portion of a solution. For problems with cable on the sea floor the last part of the solution 
may require relaxation factors on the order of 10 -3 , but may proceed quite well in the 
initial iterations with /x r ~ 10 -1 . An example of the error progress during such a problem 
is shown in figure C-l. In the upper panel, the solid line charts the error in a solution with 
a baseline relaxation of 0.2. The bottom panel shows how the relaxation factor changes as 
the solution progresses. The dashed line in the upper panel shows the error in a solution 
with a constant relaxation of 0.004 (this is the largest constant relaxation factor that 
results in a solution convergent to a = 0.001). Not only does the adaptive procedure save 
a significant amount of trial and error to determine that 0.004 is a reasonable relaxation, 
but it also reduces the number of iterations by more than a factor of two. 

C.4 Dynamic relaxation 

For certain problems, particularly moorings with cable on the bottom and very low levels 
of horizontal forcing (current and wind), the procedure outlined above can fail to converge 
when applied to the static equations. Because of the near infinite radius of curvature at the 
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Figure C-l: Error and relaxation factor during the static solution of a mooring problem 
with line on the bottom. The dashed line in the upper panel is for a solution with the 
largest constant relaxation factor (0.004) that will result in a solution convergent to an 
error level of 0.001. The solid lines are the error and relaxation factor for a solution using 
the adaptive procedure described in the text. 

touchdown point in these problems, both the shooting method initial guess (appendix D) 
and the subsequent relaxation method solution are difficult to obtain in any reasonable 
amount of computation time. A method that works well to overcome this difficulty is 
dynamic relaxation. 

In dynamic relaxation, we use the standard procedure (shooting initial guess followed 
by relaxation of the static equations) to obtain a static solution for a problem with a 
higher level of horizontal forcing. This solution is then used as the initial condition in a 
dynamic problem with the true level of current and wind applied. As time progresses in 
the dynamic simulation the mooring falls back to its true equilibrium state at the lower 
forcing level. With adaptive time-stepping, adaptive relaxation, and the physical drag 
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and damping in the problem, solutions can be obtained for significantly lower levels of 
horizontal forcing than with the standard static solution procedure. 

The procedure does not work well for three-dimensional problems in which the entire 
plane of the mooring may rotate as the true three-dimensional forces are applied. In 
these cases the time to equilibrium can be prohibitively long. Also, for either two- or 
three-dimensional cases in which the horizontal forcing is reversed from the high initial 
condition to the desired low condition, the solution can run into difficulty as the mooring 
crosses through a purely vertical configuration. 

C.5 Coordinate integration 

This solution procedure calculates the N x n dependent variables that are explicitly in¬ 
cluded in the governing equations. Because both static and dynamic governing equations 
in the formulation derived in appendix A do not explicitly include the coordinate positions 
of the nodes, these positions must be calculated in a separate procedure. The position 
of the nodes is critical to the solution; bottom boundary effects, wave forces, and current 
are all dependent on spatial position. For this reason, the coordinate positions of all the 
nodes are updated following each iteration. The integration procedure is described in 
appendix E. 
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Appendix D 


Static Initialization Procedures 


In order to solve equations A.54 through A.57 using the method outlined in appendix C 
we must calculate an initial guess at the solution. We can compute a very good first 
estimate of the solution using a shooting method to solve the governing equations without 
bending stiffness and with a simplified treatment of bottom interaction effects. Without 
bending effects the static problem reduces to two equations in two unknowns, e and cf), 
and a simple form of the shooting method can be employed. 

The shooting method solutions have the advantage that they are quite fast and provide 
a good initial solution for most problems. In many cases they are good enough to use for 
preliminary static design studies. Because of the simplifications used in these solutions, 
however, they do not provide appropriate initial conditions for the dynamic solution and 
thus we still must solve the complete static governing equations using the Newton-Raphson 
procedure described in appendix C when we want to study system dynamics. The implicit 
solutions of the complete static equations are also much more easily applied to cases 
in which the system is not singly connected - multipoint moorings and moorings with 
branches for example. 
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D.l Integration procedure 


We can derive the simplified equations from equations A.54 and A.55 by dropping the 
shear force terms, 

T'(e)— - w 0 cos <j)+ \pwdnCdt (U cos (j)+V sin (j>) | U cos <f> + V sin 0| vT T~e = 0, (D.l) 
ds 2 

Tie)— + wo sin <j)+ \p w dC dn {-U sin<£ + V cos <j))\-U sin(f) + V cos <j>\ VT+e = 0. 
os 2 

(D.2) 

Because the static boundary conditions depend on the x,y coordinates of the top and 
bottom node we also explicitly include equations for x and y. Those two equations are 

^ = (l + e)cos 4>, (D.3) 

Os 

= (1 + e) sin<£. (D.4) 

os 

For the direct integration of the governing differential equations that is inherent in the 
shooting method, this explicit inclusion of x and y is not significantly more computation¬ 
ally expensive than the integration of the coordinates after the solution that is described 
in appendix E. 

The numerical integration of the simplified governing equations proceeds from the top 
node to the anchor. Given a set of trial boundary conditions at the top node we integrate 
downwards using an explicit, fourth-order accurate Runge-Kutta algorithm [82]. If during 
the integration the calculated vertical position of a node is on or below the sea floor we 
stop the integration and assume that the remainder of the mooring is on the sea floor 
with constant angle <fi = ±f and constant tension (equal to the tension at the touchdown 
node). Sea floor slope effects are neglected in this formulation; the bottom is assumed flat 
at x — 0 . 


D.2 Iterating on the boundary conditions 

Like the outer loop iterations required to resolve the boundary conditions for the static 
solution using the relaxation procedure, shooting solutions require two levels of iteration. 
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The three types of static boundary conditions detailed in section 3.1.1 are: 

1. A buoy is on the surface, but we do not know the buoy draft. 

2. A platform is on the surface at some specified horizontal offset from the anchor. We 
do not know the mooring line tension or the angle </> at the top node. 

3. A platform is on the surface with a known tension on the mooring line. We do not 
know the angle </> at the top node. 

In all three cases we want to specify a value for strain (tension), </>, and an arbitrary 
horizontal position at the top node. We then iteratively specify a vertical position for the 
top node, integrating the simplified governing equations from node n to node 1 at each 
vertical position until the computed vertical position of the first or touchdown node is in 
fact on the bottom. The final value for the vertical coordinate of the top node is then 
used to compute corrected values for strain and <p in the outer loop. In the inner loop we 
are shooting for the correct position of the top node given some applied force. Within the 
outer loop we are shooting for the correct applied force. 

We need outer loop iterations for the relaxation procedure because the coordinate 
positions do not enter directly into the governing equations. In this case, where we do have 
coordinate positions in the governing equations, outer loop iterations are necessitated by 
the simplified treatment of the bottom. In addition to the unknown boundary conditions 
at the top node we do not know the location of the touchdown node along the mooring. 
Because the simplified bottom treatment leads to a solution for the mooring only between 
the touchdown node and the top node, the position of the touchdown point is critical. 

For the first type of boundary condition, given a guess at the draft we can calculate 
the drag and buoyancy forces and therefore T and (j) at the top node. If our first two 
guesses at the draft are the maximum and minimum available (the minimum is defined as 
the draft that will float the weight of the buoy itself and nothing more) then subsequent 
guesses can be made using bisection until the position of the top node computed in the 
inner loop corresponds to a position based on the guessed draft within some specified 
tolerance. The standard tolerance is 1% of the maximum draft. 

The third type of BC is similar in that a guess at (j) and a known T provide a complete 
force specification at the top node. With two initial trials at (j) = 0 and <p = |, the solution 
is bracketed and we can use bisection to calculate a sequence of successively better guesses 
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for (j>. The stopping procedure is the same as for the first case - the position of the top 
node computed in the inner loop must correspond with the known vertical position at the 
top within some tolerance specified as a percentage of the mesh spacing, As, at the top 
node. 

The second case is more complicated because we only know two, x and y, of the 
four boundary conditions. This requires that we iterate to find both e and <f>. Rather 
than employing a sophisticated multi-variable nonlinear root finding technique (such as 
Newton-Raphson) in the outer loop, we can use the same error correction procedure that 
we use with these boundary conditions in the outer iteration loop of the regular static 
solution. Given a guess at the vertical and horizontal components of the tension and the 
known horizontal position of the top node, we perform the inner iterations to calculate 
the vertical position of the top node. As in equation 3.7 we update the trial forces based 
on positioning error and a pseudo-“stiffness” constant, fi p , 

Fx +1 = F* - y p x k td , (D.5) 

Fy +l = Fy+ ^anchor; P-6) 


where F% y are the trial forces at iteration k in the global vertical and horizontal directions, 
respectively, x\ A is the calculated vertical coordinate of the touchdown node, and J/anchor 
is the calculated horizontal coordinate of the first node. The iterative update process is 
halted when the touchdown node is on the bottom and the anchor node is at the horizontal 
origin within a tolerance specified as a percentage of the mesh spacing at the anchor. 

The primary complication with this approach is that there is no clear best choice for the 
initial guess at the forces, F% >y , such that the iteration procedure will have a reasonable 
chance of rapid convergence. The initial forces in our implementation are based on an 
inclined catenary solution for a uniform cable with no current. Given a uniform cable 
with linear stiffness, EA, and weight per length, w 0 , the catenary solution for the position 
of the top end is 
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where F x and F y are the applied forces at the top end in the global vertical and horizontal 
directions, respectively. Given the desired position of the top end of the mooring, we use a 
two-dimensional nonlinear Newton-Raphson root finding technique to solve equations D.7 
and D.8 for F x y . In multi-segment applications, EA is calculated as the equivalent stillness 
of all segments in series, with the stiffness of each segment computed as the slope of the 
tension strain relationship at a strain of 1%. The unit weight in these cases is computed 
by summing all weight and buoyancy forces in the system and dividing by total length. 

D.3 Computing shear and curvature 

As a final step before proceeding with the relaxation solution for the complete nonlinear 
problem, the shear force, S n and curvature, G 3 are calculated numerically using centered 
differences according to equation A.56 (to calculate Q 3 ) and equation A.57 (to calculate 
S n using differences of the newly calculated Q 3 ). For boundary condition cases one and 
three, the horizontal coordinates are also translated to bring the position of the anchor to 
the origin. 


243 



244 



Appendix E 


Coordinate Integration 

Because the global coordinate variables x,y,z do not appear in any of the governing 
equations, they are integrated based on cable coordinates and cable orientation after 
each iteration of the nonlinear solver. While the coordinates do not enter directly into 
the governing equations it is important that they be updated because they are used in 
evaluating the current at a node and determining if a node is lying on the bottom. 

E.l Static solution 

For the static problem we can write differential equations for the global coordinates, x 
and y , 

(E.l) 
(E.2) 

Including these two equations directly into the static governing equations would simplify 
the handling of static boundary conditions in some cases, but only with a 50% increase 
in computational expense in the nonlinear solver. Trial implementations based on this 
approach also demonstrated convergence problems when the boundary conditions became 
part of the iterative solution. The current approach of iterating on the boundary con¬ 
ditions in a loop outside of the nonlinear solver appears to provide better stability and 
convergence. 


— = (1 + e)cos0, 

dy , . 

— = (l + e)sin <j). 
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Following the box method spatial discretization we discretize equations E.l and E.2 as 

2 (^1 — 3 i±zl\ = (1 + ej) cos 4>j + (1 + ej-i) cos (E.3) 

V s j ~ s j -1 / 

2 f——^ = (1 +ej)sm(f)j + (1 + Cj_i)sin^j_i. (E.4) 

\ s j ~ s j -1 / 

With the first node always located at the origin, we can rearrange the discretized equations 
to derive recursion relationships for xj and yj, j = 2 ... n, 

/\ • _ « 

Xj = Xj- 1 4 - j— [(1 + €j) COS <j)j + (1 + €j- 1 ) cos <f)j- 1 ] , 

A 

•_1 

Vj = Vj-i + —j— [(! + Cj) sin ^' + (1 + ej_i) sin0j_i]. 

Asj-i is the spacing between nodes j and j — 1 . 

E.2 Dynamic solution 

For the dynamic problem, we have a choice in the integration method. Equations E.l and 
E.2 are valid for the dynamic problem as are the temporal differential equations 

u l j cos 4> l j — Vj sin 4 >), (E.7) 

Uj sin 4>j + v l j cos 4>) ■ (E. 8 ) 

Either pair of equations could be incorporated into the governing equations, but again, 
only with an increase in computation expense in the nonlinear solver. With x and y 
effectively decoupled from the other six dependent variables, integration outside of the 
nonlinear solver is more efficient. 

Experience has indicated that integrating the spatial differential equations at each 
time step provides better results over long time simulations than does integration of the 
temporal equations. One explanation for this is that the spatial integration couples the 
coordinate positions of all the nodes together thus providing a strong notion of “connect¬ 
edness” at each time step. In the temporal integration the positions of the nodes are 
independent of one another, with the evolution of a node’s position in time dependent 
only on the nodal velocity and local orientation. In principle the two solutions should 


dx 

dt 

dy 

dt 


(E.5) 

(E.6) 
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be the same given the compatibility requirements enforced in the cable governing equa¬ 
tions (equations A.46 and A.47), but compatibility cannot be strictly enforced given the 
necessarily inexact solution provided by the nonlinear Newton-Raphson procedure. 

For the dynamic problem we integrate equations E.l and E.2 using the standard half¬ 
grid spatial discretization and the generalized-a method. The discretized equations are 


2 (l-a fc ) 


X S ~ x 3 ~i 


+ 2a*; 


(1 - afc) (l + e}) cos </>*• + (l + e}_j) cos 

-otk (1 + e} _1 ) cos (fir 1 + (l + e^z\) cos <f> % r\ = 0, (E.9) 


2(1 - a k ) 


yti Lil +2a jC^Ei 


- (1 - Oi k ) (l + €*•) sin (ft + (l + el.j) sin 

-Oik (1 + c} -1 ) sin^. -1 -(- (1 + e l r_\) sin <j) l r}^ = 0. (E.10) 

Rearranging terms yields the recursion relationships for the dynamic calculation of nodal 
coordinates 


x ) = x )-i + 2 (1 _' S a . fc ) ( x “ a fc) [( X + e i) cos + (! + e $-i) cos <f>U 

+ a k [(1 + e}" 1 ) cos 1 + (1 + epl) cos ^ [zf 1 - x)^ , (E.ll) 


$-i + 2(1 - ak ) t 1 “ a *) [(* + € j) sin< ^j + ( 1 + e i-i) sin ^-i] 

+ a fc [(1 + e}- 1 ) sin^}- 1 + (l + e}-l) sin $-i] “ ^ [if# -1 ~ 2/j~l] • (E.12) 
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Appendix F 


Bootstrap Monte Carlo 
Confidence Intervals 


The bootstrap method is a procedure for calculating the error in a statistically estimated 
value. The advantage to the procedure is that no assumptions are necessary about the 
underlying probability distribution of the estimated value. For the purposes of this thesis, 
the estimated values are regression coefficients calculated using least squares and the errors 
that are sought are the 95% confidence intervals of those coefficients. 

The basic procedure, as outlined by Efron and Gong [26], is as follows. Given n 
independently observed data points yi, y 2 , • • -, y n , the regression coefficients are calculated 
as 


C = / (yi,y 2 ,.-.,yn)- 


(F.l) 


C represents the best available estimate of the true value of the coefficients, C. For the 
dynamic tension model in chapter 6 y* is defined as 


yi = 




(F.2) 


The probability distribution for C is estimated using bootstrap with Monte Carlo simu¬ 
lations. A bootstrap sample, y^y^,.. • ,y*, is drawn from the original y*. The sample 
is constructed by making n random draws with replacement so that each of the original 
data points has probability 1 jn of being selected for each of the locations in the bootstrap 
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mass coefficient (kg) 

Figure F-l: m and C ( i coefficients calculated in 500 distinct bootstrap Monte Carlo sim¬ 
ulations. 95% symmetric confidence intervals for m and Cd are indicated by the dashed 
lines. The best estimates for m and Cd are indicated by the solid lines. 

sample. From this bootstrap sample we calculate a bootstrap estimate of the regression 
coefficients, 


C- = /(y;,y|,...,y;). (F.3) 

By repeating this procedure some large number of times, B, the distribution of C is 
mapped out. 

Figure F-l shows the distribution of C* = [m*, C/j for the dynamic tension model 
developed in chapter 6 with B = 500. Individual probability density functions for m and 
Cd with B = 20000 are shown in figures F-2 and F-3, respectively. With distributions for m 
and Cd there are two basic approaches to calculating confidence intervals. For equal-tailed 
confidence intervals, the bootstrap simulation results are sorted in ascending order and (for 
95% confidence) the end-points of the interval are defined as those points at indices 0.025 B 
and 0.975B in the sorted list. With this type of formulation there is equal probability 
that the true value lies above or below the interval. For mass for example, if the interval 
endpoints are defined as ml 025B and rh* 0 g75B and 8\ ~ rh — ml 02 5 £, <^2 = ^o. 975 B ~ "b 
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probability functions for m are 


< m < m + 82 ) = 0.95, 

(F.4) 

P (m < m — 81 ) = 0.025, 

(F.5) 

P (m > m + £ 2 ) = 0.025. 

(F.6) 


Hall [41] showed that symmetric intervals have coverage error O (B~ 2 ) compared to 
O (-B -1 ) for equal-tailed intervals. Because they are also slightly easier to present they 
are used throughout the thesis. The probability function for the mass coefficient with a 
symmetric confidence interval is written as 


P (m, — 8 < m < m + 8 ) = 0.95. (F.7) 

In practice, 8 is calculated using bisection. Because of the discrete nature of the count 
of points that fall in the interval, the search does not always converge to an interval 
with exactly 0.95 B points. To be conservative, the bisection algorithm always returns an 
interval that contains at least 0.951? points. 


251 




Figure F-2: Probability density function for m based on 20000 distinct bootstrap Monte 
Carlo simulations. The best estimate for m from the original least squares fit is indi¬ 
cated by the solid vertical line. The mean from the bootstrap simulations is indicated by 
the dashed vertical line. 95% symmetric [ ] and equal-tailed ( ) confidence intervals are 
indicated on the horizontal axis. 



Figure F-3: Probability density function for Cd based on 20000 distinct bootstrap Monte 
Carlo simulations. Other markings are the same as in figure F-2. 
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Appendix G 


Catenary formulae 


For an inextensible line with no current, and vertical and horizontal forces applied at the 
top point, F v and Fh, respectively, the catenary expression for the vertical coordinate of 
the top point is 


z — 


Wq 



(G.l) 


With excess scope remaining on the bottom the vertical force at the top of the mooring, 
F v , must equal the suspended weight, w 0 L. The total tension, T, at the top is simply 
\Jf% + F 2 . If the top of the mooring is at the surface, z = H, depth, tension and length 
are related by 


T = woH+s]t*-{ W oL) 2 . 


(G.2) 


After some manipulation, the non-dimensional mean tension can be written as 

1 

T ~ 2 

Ar = - 


, L 2 \ 



(G.3) 

T 2 \ 


H 2 1 ) ‘ 

(G.4) 


The suspended length as a function of non-dimensional mean tension is 


L — HV2At + 1. 


(G.5) 
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The rate of increase of the suspended length with increasing At is 

dr L 

This rate slows as the scope of the mooring increases. 


(G.6) 
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